From 83dd99541ffbc2434d3cc71300c646b02e2f46d5 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Thu, 12 Dec 2024 14:54:06 -0800 Subject: [PATCH 01/17] begin composite and size distribution --- src/pybamm/expression_tree/averages.py | 15 ++- .../full_battery_models/base_battery_model.py | 1 - .../submodels/interface/base_interface.py | 29 ++++-- .../interface/kinetics/base_kinetics.py | 12 ++- .../interface/lithium_plating/base_plating.py | 21 ++++ .../interface/lithium_plating/no_plating.py | 25 ++++- .../open_circuit_potential/base_ocp.py | 9 +- .../current_sigmoid_ocp.py | 2 +- .../open_circuit_potential/single_ocp.py | 4 +- .../submodels/interface/sei/base_sei.py | 18 +++- .../models/submodels/interface/sei/no_sei.py | 10 ++ .../size_distribution_parameters.py | 95 ++++++++++++++----- 12 files changed, 192 insertions(+), 49 deletions(-) diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py index 11538ea153..75c68392c5 100644 --- a/src/pybamm/expression_tree/averages.py +++ b/src/pybamm/expression_tree/averages.py @@ -306,6 +306,10 @@ def xyz_average(symbol: pybamm.Symbol) -> pybamm.Symbol: return yz_average(x_average(symbol)) +def xyzs_average(symbol: pybamm.Symbol) -> pybamm.Symbol: + return xyz_average(size_average(symbol)) + + def r_average(symbol: pybamm.Symbol) -> pybamm.Symbol: """ Convenience function for creating an average in the r-direction. @@ -373,9 +377,18 @@ def size_average( # If symbol doesn't have a domain, or doesn't have "negative particle size" # or "positive particle size" as a domain, it's average value is itself if symbol.domain == [] or not any( - domain in [["negative particle size"], ["positive particle size"]] + domain + in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ] for domain in list(symbol.domains.values()) ): + print(f"symbol: {symbol.domain}") return symbol # If symbol is a primary broadcast to "particle size", take the orphan diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py index fd404cec75..6fa110545e 100644 --- a/src/pybamm/models/full_battery_models/base_battery_model.py +++ b/src/pybamm/models/full_battery_models/base_battery_model.py @@ -623,7 +623,6 @@ def __init__(self, extra_options): if options["particle phases"] not in ["1", ("1", "1")]: if not ( options["surface form"] != "false" - and options["particle size"] == "single" and options["particle"] == "Fickian diffusion" ): raise pybamm.OptionError( diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index 0ad08d5454..8978343bba 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -46,6 +46,8 @@ def __init__(self, param, domain, reaction, options, phase="primary"): self.reaction_name = self.phase_name + self.reaction_name self.reaction = reaction + domain_options = getattr(self.options, domain) + self.size_distribution = domain_options["particle size"] == "distribution" def _get_exchange_current_density(self, variables): """ @@ -93,8 +95,10 @@ def _get_exchange_current_density(self, variables): # "current collector" c_e = pybamm.PrimaryBroadcast(c_e, ["current collector"]) # broadcast c_e, T onto "particle size" - c_e = pybamm.PrimaryBroadcast(c_e, [f"{domain} particle size"]) - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + c_e = pybamm.PrimaryBroadcast( + c_e, [f"{domain} {phase_name}particle size"] + ) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: c_s_surf = variables[ @@ -215,7 +219,14 @@ def _get_standard_interfacial_current_variables(self, j): # Size average. For j variables that depend on particle size, see # "_get_standard_size_distribution_interfacial_current_variables" - if j.domain in [["negative particle size"], ["positive particle size"]]: + if j.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: j = pybamm.size_average(j) # Average, and broadcast if necessary j_av = pybamm.x_average(j) @@ -299,9 +310,15 @@ def _get_standard_volumetric_current_density_variables(self, variables): f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]" ] - j = variables[ - f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]" - ] + if False: + j = variables[ + f"{Domain} electrode {phase_name}{reaction_name}interfacial current density distribution [A.m-2]" + ] + else: + j = variables[ + f"{Domain} electrode {phase_name}interfacial current density [A.m-2]" + ] + print(reaction_name) a_j = a * j a_j_av = pybamm.x_average(a_j) diff --git a/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py index 7cfa83631b..c3b6ed5329 100644 --- a/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py +++ b/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py @@ -76,7 +76,9 @@ def get_coupled_variables(self, variables): self.reaction == "lithium-ion main" and domain_options["particle size"] == "distribution" ): - delta_phi = pybamm.PrimaryBroadcast(delta_phi, [f"{domain} particle size"]) + delta_phi = pybamm.PrimaryBroadcast( + delta_phi, [f"{domain} {phase_name}particle size"] + ) # Get exchange-current density. For MSMR models we calculate the exchange # current density for each reaction, then sum these to give a total exchange @@ -169,7 +171,7 @@ def get_coupled_variables(self, variables): elif j0.domain in ["current collector", ["current collector"]]: T = variables["X-averaged cell temperature [K]"] u = variables[f"X-averaged {domain} electrode interface utilisation"] - elif j0.domain == [f"{domain} particle size"]: + elif j0.domain == [f"{domain} {phase_name}particle size"]: if j0.domains["secondary"] != [f"{domain} electrode"]: T = variables["X-averaged cell temperature [K]"] u = variables[f"X-averaged {domain} electrode interface utilisation"] @@ -178,7 +180,7 @@ def get_coupled_variables(self, variables): u = variables[f"{Domain} electrode interface utilisation"] # Broadcast T onto "particle size" domain - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: T = variables[f"{Domain} electrode temperature [K]"] u = variables[f"{Domain} electrode interface utilisation"] @@ -198,7 +200,9 @@ def get_coupled_variables(self, variables): else: j = self._get_kinetics(j0, ne, eta_r, T, u) - if j.domain == [f"{domain} particle size"]: + if j.domain == [f"{domain} particle size"] or j.domain == [ + f"{domain} {phase_name}particle size" + ]: # If j depends on particle size, get size-dependent "distribution" # variables first variables.update( diff --git a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py index 0530476f9d..e9a2b64b22 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py @@ -137,3 +137,24 @@ def _get_standard_reaction_variables(self, j_stripping): } return variables + + def _get_standard_size_distribution_reaction_variables(self, j_stripping): + """ + A private function to obtain the standard variables which + can be derived from the lithum stripping interfacial reaction current + """ + domain, Domain = self.domain_Domain + phase_name = self.phase_name + j_stripping_sav = pybamm.size_average(j_stripping) + j_stripping_av = pybamm.x_average(j_stripping_sav) + + variables = { + f"{Domain} electrode {phase_name}lithium plating " + "interfacial current density distribution [A.m-2]": j_stripping, + f"{Domain} electrode {phase_name}lithium plating " + "interfacial current density [A.m-2]": j_stripping_sav, + f"X-averaged {domain} electrode {phase_name}lithium plating " + "interfacial current density [A.m-2]": j_stripping_av, + } + + return variables diff --git a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py index ddc13cfb97..00bd0bf164 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py @@ -20,12 +20,27 @@ def __init__(self, param, domain, options=None, phase="primary"): super().__init__(param, domain, options=options, phase=phase) def get_fundamental_variables(self): - zero = pybamm.FullBroadcast( - pybamm.Scalar(0), f"{self.domain} electrode", "current collector" - ) + phase_name = self.phase_name + if self.size_distribution is False: + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), f"{self.domain} electrode", "current collector" + ) + else: + zero = pybamm.PrimaryBroadcast( + pybamm.Scalar(0), f"{self.domain} {phase_name}particle size" + ) variables = self._get_standard_concentration_variables(zero, zero) - variables.update(self._get_standard_overpotential_variables(zero)) - variables.update(self._get_standard_reaction_variables(zero)) + if self.size_distribution is False: + variables.update(self._get_standard_overpotential_variables(zero)) + variables.update(self._get_standard_reaction_variables(zero)) + else: + variables.update( + self._get_standard_size_distribution_overpotential_variables(zero) + ) + variables.update( + self._get_standard_size_distribution_reaction_variables(zero) + ) + return variables def get_coupled_variables(self, variables): diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py index 309eb343c0..e2c4a43551 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py @@ -49,7 +49,14 @@ def _get_standard_ocp_variables(self, ocp_surf, ocp_bulk, dUdT): reaction_name = self.reaction_name # Update size variables then size average. - if ocp_surf.domain in [["negative particle size"], ["positive particle size"]]: + if ocp_surf.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: variables = self._get_standard_size_distribution_ocp_variables( ocp_surf, dUdT ) diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py index fa55ba1079..8fc41d1711 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py @@ -36,7 +36,7 @@ def get_coupled_variables(self, variables): ): sto_surf = sto_surf.orphans[0] T = T.orphans[0] - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: sto_surf = variables[ f"{Domain} {phase_name}particle surface stoichiometry" diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py index 932e2e3eff..42a4b53ebb 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py @@ -25,7 +25,7 @@ def get_coupled_variables(self, variables): ): sto_surf = sto_surf.orphans[0] T = T.orphans[0] - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: sto_surf = variables[ f"{Domain} {phase_name}particle surface stoichiometry" @@ -42,7 +42,7 @@ def get_coupled_variables(self, variables): # Bulk OCP is from the average SOC and temperature sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"] - T_bulk = pybamm.xyz_average(pybamm.size_average(T)) + T_bulk = pybamm.xyzs_average(T) ocp_bulk = self.phase_param.U(sto_bulk, T_bulk) elif self.reaction == "lithium metal plating": T = variables[f"{Domain} electrode temperature [K]"] diff --git a/src/pybamm/models/submodels/interface/sei/base_sei.py b/src/pybamm/models/submodels/interface/sei/base_sei.py index 8a4176d416..ee8d6df7e3 100644 --- a/src/pybamm/models/submodels/interface/sei/base_sei.py +++ b/src/pybamm/models/submodels/interface/sei/base_sei.py @@ -152,8 +152,13 @@ def _get_standard_concentration_variables(self, variables): L = variables[f"{Domain} {reaction_name}thickness [m]"] n_SEI = L * L_to_n - n_SEI_xav = pybamm.x_average(n_SEI) - n_SEI_av = pybamm.yz_average(n_SEI_xav) + if self.size_distribution: + n_SEI_sav = pybamm.size_average(n_SEI) + n_SEI_xav = pybamm.x_average(n_SEI_sav) + n_SEI_av = pybamm.yz_average(n_SEI_xav) + else: + n_SEI_xav = pybamm.x_average(n_SEI) + n_SEI_av = pybamm.yz_average(n_SEI_xav) # Calculate change in SEI concentration with respect to initial state delta_n_SEI = n_SEI_av - n_SEI_0 @@ -188,8 +193,13 @@ def _get_standard_concentration_variables(self, variables): roughness = variables[f"{Domain} electrode roughness ratio"] n_SEI_cr = L_cr * L_to_n * (roughness - 1) # SEI on cracks concentration - n_SEI_cr_xav = pybamm.x_average(n_SEI_cr) - n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) + if self.size_distribution: + n_SEI_cr_sav = pybamm.size_average(n_SEI_cr) + n_SEI_cr_xav = pybamm.x_average(n_SEI_cr_sav) + n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) + else: + n_SEI_cr_xav = pybamm.x_average(n_SEI_cr) + n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) # Calculate change in SEI cracks concentration # Initial state depends on roughness (to avoid division by zero) diff --git a/src/pybamm/models/submodels/interface/sei/no_sei.py b/src/pybamm/models/submodels/interface/sei/no_sei.py index 27cea3b69b..242c8676ea 100644 --- a/src/pybamm/models/submodels/interface/sei/no_sei.py +++ b/src/pybamm/models/submodels/interface/sei/no_sei.py @@ -32,12 +32,22 @@ def get_fundamental_variables(self): domain = self.domain.lower() if self.reaction_loc == "interface": zero = pybamm.PrimaryBroadcast(pybamm.Scalar(0), "current collector") + elif self.size_distribution: + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), + f"{domain} {self.phase_name}particle size", + {"secondary": f"{domain} electrode", "tertiary": "current collector"}, + ) else: zero = pybamm.FullBroadcast( pybamm.Scalar(0), f"{domain} electrode", "current collector" ) variables = self._get_standard_thickness_variables(zero) variables.update(self._get_standard_reaction_variables(zero)) + if self.size_distribution: + variables.update( + self._get_standard_size_distribution_interfacial_current_variables(zero) + ) return variables def get_coupled_variables(self, variables): diff --git a/src/pybamm/parameters/size_distribution_parameters.py b/src/pybamm/parameters/size_distribution_parameters.py index 60383fff50..6693cd0af3 100644 --- a/src/pybamm/parameters/size_distribution_parameters.py +++ b/src/pybamm/parameters/size_distribution_parameters.py @@ -9,14 +9,23 @@ def get_size_distribution_parameters( param, R_n_av=None, + R_n_av_prim=None, + R_n_av_sec=None, R_p_av=None, sd_n=0.3, + sd_n_prim=0.3, + sd_n_sec=0.3, sd_p=0.3, R_min_n=None, + R_min_n_prim=None, + R_min_n_sec=None, R_min_p=None, R_max_n=None, + R_max_n_prim=None, + R_max_n_sec=None, R_max_p=None, working_electrode="both", + composite="neither", ): """ A convenience method to add standard area-weighted particle-size distribution @@ -62,30 +71,68 @@ def get_size_distribution_parameters( """ if working_electrode == "both": # Radii from given parameter set - R_n_typ = param["Negative particle radius [m]"] - - # Set the mean particle radii - R_n_av = R_n_av or R_n_typ - - # Minimum radii - R_min_n = R_min_n or np.max([0, 1 - sd_n * 5]) - - # Max radii - R_max_n = R_max_n or (1 + sd_n * 5) - - # Area-weighted particle-size distribution - def f_a_dist_n(R): - return lognormal(R, R_n_av, sd_n * R_n_av) - - param.update( - { - "Negative minimum particle radius [m]": R_min_n * R_n_av, - "Negative maximum particle radius [m]": R_max_n * R_n_av, - "Negative area-weighted " - + "particle-size distribution [m-1]": f_a_dist_n, - }, - check_already_exists=False, - ) + if composite != "negative" and composite != "both": + R_n_typ = param["Negative particle radius [m]"] + # Set the mean particle radii + R_n_av = R_n_av or R_n_typ + # Minimum radii + R_min_n = R_min_n or np.max([0, 1 - sd_n * 5]) + + # Max radii + R_max_n = R_max_n or (1 + sd_n * 5) + + # Area-weighted particle-size distribution + def f_a_dist_n(R): + return lognormal(R, R_n_av, sd_n * R_n_av) + + param.update( + { + "Negative minimum particle radius [m]": R_min_n * R_n_av, + "Negative maximum particle radius [m]": R_max_n * R_n_av, + "Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n, + }, + check_already_exists=False, + ) + else: + R_n_typ_prim = param["Primary: Negative particle radius [m]"] + R_n_typ_sec = param["Secondary: Negative particle radius [m]"] + # Set the mean particle radii + R_n_av_prim = R_n_av_prim or R_n_typ_prim + R_n_av_sec = R_n_av_sec or R_n_typ_sec + # Minimum radii + R_min_n_prim = R_min_n_prim or np.max([0, 1 - sd_n_prim * 5]) + R_min_n_sec = R_min_n_sec or np.max([0, 1 - sd_n_sec * 5]) + + # Max radii + R_max_n_prim = R_max_n_prim or (1 + sd_n_prim * 5) + R_max_n_sec = R_max_n_sec or (1 + sd_n_sec * 5) + + # Area-weighted particle-size distribution + def f_a_dist_n_prim(R): + return lognormal(R, R_n_av_prim, sd_n_prim * R_n_av_prim) + + def f_a_dist_n_sec(R): + return lognormal(R, R_n_av_sec, sd_n_sec * R_n_av_sec) + + param.update( + { + "Primary: Negative minimum particle radius [m]": R_min_n_prim + * R_n_av_prim, + "Primary: Negative maximum particle radius [m]": R_max_n_prim + * R_n_av_prim, + "Primary: Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n_prim, + "Secondary: Negative minimum particle radius [m]": R_min_n_sec + * R_n_av_sec, + "Secondary: Negative maximum particle radius [m]": R_max_n_sec + * R_n_av_sec, + "Secondary: Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n_sec, + }, + check_already_exists=False, + ) + # Radii from given parameter set R_p_typ = param["Positive particle radius [m]"] From 843a9e160c89b30681922cf81b0e1e7c9a2503e4 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Fri, 13 Dec 2024 11:26:58 -0800 Subject: [PATCH 02/17] add phase prefactor; still not working --- src/pybamm/expression_tree/averages.py | 12 +++++++++-- .../submodels/particle/base_particle.py | 21 ++++++++++--------- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py index 75c68392c5..4866ad78c7 100644 --- a/src/pybamm/expression_tree/averages.py +++ b/src/pybamm/expression_tree/averages.py @@ -410,10 +410,18 @@ def size_average( R = pybamm.SpatialVariable( "R", domains=symbol.domains, coord_sys="cartesian" ) - if ["negative particle size"] in symbol.domains.values(): + if ["negative particle size"] in symbol.domains.values() or [ + "negative primary particle size" + ] in symbol.domains.values(): f_a_dist = geo.n.prim.f_a_dist(R) - elif ["positive particle size"] in symbol.domains.values(): + elif ["negative secondary particle size"] in symbol.domains.values(): + f_a_dist = geo.n.sec.f_a_dist(R) + elif ["positive particle size"] in symbol.domains.values() or [ + "positive primary particle size" + ] in symbol.domains.values(): f_a_dist = geo.p.prim.f_a_dist(R) + elif ["positive secondary particle size"] in symbol.domains.values(): + f_a_dist = geo.p.sec.f_a_dist(R) return SizeAverage(symbol, f_a_dist) diff --git a/src/pybamm/models/submodels/particle/base_particle.py b/src/pybamm/models/submodels/particle/base_particle.py index b774e58a0c..b5b8fb4499 100644 --- a/src/pybamm/models/submodels/particle/base_particle.py +++ b/src/pybamm/models/submodels/particle/base_particle.py @@ -185,6 +185,7 @@ def _get_distribution_variables(self, R): """ domain, Domain = self.domain_Domain phase_name = self.phase_name + Phase_prefactor = self.phase_param.phase_prefactor R_typ = self.phase_param.R_typ # [m] # Particle-size distribution (area-weighted) @@ -237,21 +238,21 @@ def _get_distribution_variables(self, R): variables = { f"{Domain} {phase_name}particle sizes": R / R_typ, - f"{Domain} {phase_name}particle sizes [m]": R, - f"{Domain} area-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} {phase_name}particle sizes [m]": R, + f"{Phase_prefactor}{Domain} area-weighted {phase_name}particle-size" " distribution [m-1]": f_a_dist, - f"{Domain} volume-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}particle-size" " distribution [m-1]": f_v_dist, - f"{Domain} number-based {phase_name}particle-size" + f"{Phase_prefactor}{Domain} number-based {phase_name}particle-size" " distribution [m-1]": f_num_dist, - f"{Domain} area-weighted mean particle radius [m]": R_a_mean, - f"{Domain} volume-weighted mean particle radius [m]": R_v_mean, - f"{Domain} number-based mean particle radius [m]": R_num_mean, - f"{Domain} area-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} area-weighted mean particle radius [m]": R_a_mean, + f"{Phase_prefactor}{Domain} volume-weighted mean particle radius [m]": R_v_mean, + f"{Phase_prefactor}{Domain} number-based mean particle radius [m]": R_num_mean, + f"{Phase_prefactor}{Domain} area-weighted {phase_name}particle-size" " standard deviation [m]": sd_a, - f"{Domain} volume-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}particle-size" " standard deviation [m]": sd_v, - f"{Domain} number-based {phase_name}particle-size" + f"{Phase_prefactor}{Domain} number-based {phase_name}particle-size" " standard deviation [m]": sd_num, # X-averaged sizes and distributions f"X-averaged {domain} {phase_name}particle sizes [m]": pybamm.x_average(R), From bc1a2b8894be635426cda8b3ad7943d78db52900 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Tue, 17 Dec 2024 11:24:35 -0800 Subject: [PATCH 03/17] particle size and composite electrode works --- src/pybamm/expression_tree/averages.py | 29 ++++++++++++-- src/pybamm/geometry/battery_geometry.py | 24 ++++++++++++ src/pybamm/geometry/standard_spatial_vars.py | 39 +++++++++++++++++++ src/pybamm/meshes/meshes.py | 1 + .../full_battery_models/base_battery_model.py | 12 ++++++ .../active_material/base_active_material.py | 14 ++++++- .../interface/lithium_plating/base_plating.py | 8 +++- .../interface/lithium_plating/no_plating.py | 9 ++++- .../submodels/particle/fickian_diffusion.py | 7 ++-- src/pybamm/parameters/base_parameters.py | 4 ++ src/pybamm/parameters/geometric_parameters.py | 6 +++ src/pybamm/parameters/parameter_values.py | 1 + 12 files changed, 143 insertions(+), 11 deletions(-) diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py index 4866ad78c7..33f38c68b3 100644 --- a/src/pybamm/expression_tree/averages.py +++ b/src/pybamm/expression_tree/averages.py @@ -407,9 +407,32 @@ def size_average( else: if f_a_dist is None: geo = pybamm.geometric_parameters - R = pybamm.SpatialVariable( - "R", domains=symbol.domains, coord_sys="cartesian" - ) + if "negative" in symbol.domain[0]: + if "primary" in symbol.domain[0]: + R = pybamm.SpatialVariable( + "R_n_prim", domains=symbol.domains, coord_sys="cartesian" + ) + elif "secondary" in symbol.domain[0]: + R = pybamm.SpatialVariable( + "R_n_sec", domains=symbol.domains, coord_sys="cartesian" + ) + else: + R = pybamm.SpatialVariable( + "R_n", domains=symbol.domains, coord_sys="cartesian" + ) + else: + if "primary" in symbol.domain[0]: + R = pybamm.SpatialVariable( + "R_p_prim", domains=symbol.domains, coord_sys="cartesian" + ) + elif "secondary" in symbol.domain[0]: + R = pybamm.SpatialVariable( + "R_p_sec", domains=symbol.domains, coord_sys="cartesian" + ) + else: + R = pybamm.SpatialVariable( + "R_p", domains=symbol.domains, coord_sys="cartesian" + ) if ["negative particle size"] in symbol.domains.values() or [ "negative primary particle size" ] in symbol.domains.values(): diff --git a/src/pybamm/geometry/battery_geometry.py b/src/pybamm/geometry/battery_geometry.py index 9be08ff619..3e1986507a 100644 --- a/src/pybamm/geometry/battery_geometry.py +++ b/src/pybamm/geometry/battery_geometry.py @@ -83,6 +83,18 @@ def battery_geometry( "negative particle size": {"R_n": {"min": R_min_n, "max": R_max_n}}, } ) + phases = int(options.negative["particle phases"]) + if phases >= 2: + geometry.update( + { + "negative primary particle size": { + "R_n_prim": {"min": R_min_n, "max": R_max_n} + }, + "negative secondary particle size": { + "R_n_sec": {"min": R_min_n, "max": R_max_n} + }, + } + ) if ( options is not None and options.positive["particle size"] == "distribution" @@ -95,6 +107,18 @@ def battery_geometry( "positive particle size": {"R_p": {"min": R_min_p, "max": R_max_p}}, } ) + phases = int(options.positive["particle phases"]) + if phases >= 2: + geometry.update( + { + "positive primary particle size": { + "R_p_prim": {"min": R_min_p, "max": R_max_p} + }, + "positive secondary particle size": { + "R_p_sec": {"min": R_min_p, "max": R_max_p} + }, + } + ) # Add current collector domains current_collector_dimension = options["dimensionality"] diff --git a/src/pybamm/geometry/standard_spatial_vars.py b/src/pybamm/geometry/standard_spatial_vars.py index 565908deb9..b6638694f7 100644 --- a/src/pybamm/geometry/standard_spatial_vars.py +++ b/src/pybamm/geometry/standard_spatial_vars.py @@ -98,6 +98,7 @@ }, coord_sys="cartesian", ) + R_p = pybamm.SpatialVariable( "R_p", domain=["positive particle size"], @@ -108,6 +109,44 @@ coord_sys="cartesian", ) +R_n_prim = pybamm.SpatialVariable( + "R_n_prim", + domain=["negative primary particle size"], + auxiliary_domains={ + "secondary": "negative electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) +R_p_prim = pybamm.SpatialVariable( + "R_p_prim", + domain=["positive primary particle size"], + auxiliary_domains={ + "secondary": "positive electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) + +R_n_sec = pybamm.SpatialVariable( + "R_n_sec", + domain=["negative secondary particle size"], + auxiliary_domains={ + "secondary": "negative electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) +R_p_sec = pybamm.SpatialVariable( + "R_p_sec", + domain=["positive secondary particle size"], + auxiliary_domains={ + "secondary": "positive electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) + # Domains at cell edges x_n_edge = pybamm.SpatialVariableEdge( "x_n", diff --git a/src/pybamm/meshes/meshes.py b/src/pybamm/meshes/meshes.py index 3ec291b1b8..45006bf88a 100644 --- a/src/pybamm/meshes/meshes.py +++ b/src/pybamm/meshes/meshes.py @@ -178,6 +178,7 @@ def combine_submeshes(self, *submeshnames): raise pybamm.DomainError( "trying to combine two meshes in different coordinate systems" ) + print(submeshnames) combined_submesh_edges = np.concatenate( [self[submeshnames[0]].edges] + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]] diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py index 6fa110545e..e47a6c05c8 100644 --- a/src/pybamm/models/full_battery_models/base_battery_model.py +++ b/src/pybamm/models/full_battery_models/base_battery_model.py @@ -867,6 +867,10 @@ def default_var_pts(self): "z": 10, "R_n": 30, "R_p": 30, + "R_n_prim": 30, + "R_p_prim": 30, + "R_n_sec": 30, + "R_p_sec": 30, } # Reduce the default points for 2D current collectors if self.options["dimensionality"] == 2: @@ -887,6 +891,10 @@ def default_submesh_types(self): "positive secondary particle": pybamm.Uniform1DSubMesh, "negative particle size": pybamm.Uniform1DSubMesh, "positive particle size": pybamm.Uniform1DSubMesh, + "negative primary particle size": pybamm.Uniform1DSubMesh, + "positive primary particle size": pybamm.Uniform1DSubMesh, + "negative secondary particle size": pybamm.Uniform1DSubMesh, + "positive secondary particle size": pybamm.Uniform1DSubMesh, } if self.options["dimensionality"] == 0: base_submeshes["current collector"] = pybamm.SubMesh0D @@ -909,6 +917,10 @@ def default_spatial_methods(self): "positive secondary particle": pybamm.FiniteVolume(), "negative particle size": pybamm.FiniteVolume(), "positive particle size": pybamm.FiniteVolume(), + "negative primary particle size": pybamm.FiniteVolume(), + "positive primary particle size": pybamm.FiniteVolume(), + "negative secondary particle size": pybamm.FiniteVolume(), + "positive secondary particle size": pybamm.FiniteVolume(), } if self.options["dimensionality"] == 0: # 0D submesh - use base spatial method diff --git a/src/pybamm/models/submodels/active_material/base_active_material.py b/src/pybamm/models/submodels/active_material/base_active_material.py index ea0e826e09..4827c1f6d4 100644 --- a/src/pybamm/models/submodels/active_material/base_active_material.py +++ b/src/pybamm/models/submodels/active_material/base_active_material.py @@ -81,9 +81,19 @@ def _get_standard_active_material_variables(self, eps_solid): R = self.phase_param.R elif domain_options["particle size"] == "distribution": if self.domain == "negative": - R_ = pybamm.standard_spatial_vars.R_n + if self.phase_name == "primary ": + R_ = pybamm.standard_spatial_vars.R_n_prim + elif self.phase_name == "secondary ": + R_ = pybamm.standard_spatial_vars.R_n_sec + else: + R_ = pybamm.standard_spatial_vars.R_n elif self.domain == "positive": - R_ = pybamm.standard_spatial_vars.R_p + if self.phase_name == "primary ": + R_ = pybamm.standard_spatial_vars.R_p_prim + elif self.phase_name == "secondary ": + R_ = pybamm.standard_spatial_vars.R_p_sec + else: + R_ = pybamm.standard_spatial_vars.R_p R = pybamm.size_average(R_) R_av = pybamm.x_average(R) diff --git a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py index e9a2b64b22..9425e375b2 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py @@ -64,7 +64,9 @@ def _get_standard_concentration_variables(self, c_plated_Li, c_dead_Li): phase_name = self.phase_name phase_param = self.phase_param domain, Domain = self.domain_Domain - + if self.size_distribution is True: + c_plated_Li = pybamm.size_average(c_plated_Li) + c_dead_Li = pybamm.size_average(c_dead_Li) # Set scales to one for the "no plating" model so that they are not required # by parameter values in general if isinstance(self, pybamm.lithium_plating.NoPlating): @@ -88,6 +90,10 @@ def _get_standard_concentration_variables(self, c_plated_Li, c_dead_Li): L_dead_Li_av = pybamm.x_average(L_dead_Li) Q_dead_Li = c_dead_Li_av * L_k * self.param.L_y * self.param.L_z + # if self.size_distribution is True: + # Q_plated_Li = pybamm.size_average(Q_plated_Li) + # Q_dead_Li = pybamm.size_average(Q_dead_Li) + variables = { f"{Domain} {phase_name}lithium plating concentration " "[mol.m-3]": c_plated_Li, diff --git a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py index 00bd0bf164..64af5876df 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py @@ -26,8 +26,13 @@ def get_fundamental_variables(self): pybamm.Scalar(0), f"{self.domain} electrode", "current collector" ) else: - zero = pybamm.PrimaryBroadcast( - pybamm.Scalar(0), f"{self.domain} {phase_name}particle size" + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), + f"{self.domain} {phase_name}particle size", + { + "secondary": f"{self.domain} electrode", + "tertiary": "current collector", + }, ) variables = self._get_standard_concentration_variables(zero, zero) if self.size_distribution is False: diff --git a/src/pybamm/models/submodels/particle/fickian_diffusion.py b/src/pybamm/models/submodels/particle/fickian_diffusion.py index 634e2ce730..5ae5d3a4da 100644 --- a/src/pybamm/models/submodels/particle/fickian_diffusion.py +++ b/src/pybamm/models/submodels/particle/fickian_diffusion.py @@ -31,7 +31,7 @@ def __init__(self, param, domain, options, phase="primary", x_average=False): def get_fundamental_variables(self): domain, Domain = self.domain_Domain phase_name = self.phase_name - + Phase_prefactor = self.phase_param.phase_prefactor variables = {} if self.size_distribution is False: if self.x_average is False: @@ -81,7 +81,7 @@ def get_fundamental_variables(self): ) variables = self._get_distribution_variables(R) f_v_dist = variables[ - f"{Domain} volume-weighted {phase_name}" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}" "particle-size distribution [m-1]" ] else: @@ -130,6 +130,7 @@ def get_fundamental_variables(self): def get_coupled_variables(self, variables): domain, Domain = self.domain_Domain phase_name = self.phase_name + Phase_prefactor = self.phase_param.phase_prefactor if self.size_distribution is False: if self.x_average is False: @@ -219,7 +220,7 @@ def get_coupled_variables(self, variables): ) variables.update(self._get_standard_flux_distribution_variables(N_s)) # Size-averaged flux variables - R = variables[f"{Domain} {phase_name}particle sizes [m]"] + R = variables[f"{Phase_prefactor}{Domain} {phase_name}particle sizes [m]"] f_a_dist = self.phase_param.f_a_dist(R) D_eff = pybamm.Integral(f_a_dist * D_eff, R) N_s = pybamm.Integral(f_a_dist * N_s, R) diff --git a/src/pybamm/parameters/base_parameters.py b/src/pybamm/parameters/base_parameters.py index c686665019..246e7e5af7 100644 --- a/src/pybamm/parameters/base_parameters.py +++ b/src/pybamm/parameters/base_parameters.py @@ -88,6 +88,8 @@ def set_phase_name(self): self.phase == "primary" and getattr(self.main_param.options, self.domain)["particle phases"] == "1" ): + print("hello old chap") + print(self.domain_Domain) # Only one phase, no need to distinguish between # "primary" and "secondary" self.phase_name = "" @@ -95,6 +97,8 @@ def set_phase_name(self): else: # add a space so that we can use "" or (e.g.) "primary " interchangeably # when naming variables + print("hello young fellow") + print(self.domain_Domain, self.phase) self.phase_name = self.phase + " " self.phase_prefactor = self.phase.capitalize() + ": " diff --git a/src/pybamm/parameters/geometric_parameters.py b/src/pybamm/parameters/geometric_parameters.py index 0155bbe03c..415ffb8e72 100644 --- a/src/pybamm/parameters/geometric_parameters.py +++ b/src/pybamm/parameters/geometric_parameters.py @@ -165,6 +165,12 @@ def f_a_dist(self, R): Dimensional electrode area-weighted particle-size distribution """ Domain = self.domain.capitalize() + if "primary" in R.domains["primary"][0]: + self.phase_prefactor = "Primary: " + elif "secondary" in R.domains["primary"][0]: + self.phase_prefactor = "Secondary: " + else: + pass inputs = {f"{self.phase_prefactor}{Domain} particle-size variable [m]": R} return pybamm.FunctionParameter( f"{self.phase_prefactor}{Domain} " diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py index 9ddbfd2fd3..318184c356 100644 --- a/src/pybamm/parameters/parameter_values.py +++ b/src/pybamm/parameters/parameter_values.py @@ -666,6 +666,7 @@ def _process_symbol(self, symbol): elif isinstance(symbol, pybamm.FunctionParameter): function_name = self[symbol.name] + print(symbol.name) if isinstance( function_name, (numbers.Number, pybamm.Interpolant, pybamm.InputParameter), From d1cc2f12f8d779a7290c3e8b07bd7367f72bf1f0 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Tue, 17 Dec 2024 11:36:30 -0800 Subject: [PATCH 04/17] remove debug print statement --- src/pybamm/parameters/base_parameters.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/pybamm/parameters/base_parameters.py b/src/pybamm/parameters/base_parameters.py index 246e7e5af7..c686665019 100644 --- a/src/pybamm/parameters/base_parameters.py +++ b/src/pybamm/parameters/base_parameters.py @@ -88,8 +88,6 @@ def set_phase_name(self): self.phase == "primary" and getattr(self.main_param.options, self.domain)["particle phases"] == "1" ): - print("hello old chap") - print(self.domain_Domain) # Only one phase, no need to distinguish between # "primary" and "secondary" self.phase_name = "" @@ -97,8 +95,6 @@ def set_phase_name(self): else: # add a space so that we can use "" or (e.g.) "primary " interchangeably # when naming variables - print("hello young fellow") - print(self.domain_Domain, self.phase) self.phase_name = self.phase + " " self.phase_prefactor = self.phase.capitalize() + ": " From 6b303a87f05cabfac031e22dcb3d31f98eacbde4 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Tue, 17 Dec 2024 12:54:31 -0800 Subject: [PATCH 05/17] fix something --- src/pybamm/models/submodels/interface/base_interface.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index 8978343bba..c511d4ded0 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -427,9 +427,11 @@ def _get_standard_size_distribution_interfacial_current_variables(self, j): if j.domains["secondary"] == [f"{domain} electrode"]: # x-average j_xav = pybamm.x_average(j) - else: + elif getattr(self, "reaction_loc", None) != "interface": j_xav = j j = pybamm.SecondaryBroadcast(j_xav, [f"{domain} electrode"]) + else: + j_xav = j variables = { f"{Domain} electrode {reaction_name}" From 5ca70ae24120c2726e26089b91c2fb2b4c67d60a Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Tue, 17 Dec 2024 16:21:53 -0800 Subject: [PATCH 06/17] almost working; half cell still doesnt work --- .../submodels/interface/base_interface.py | 16 ++++++--- .../submodels/interface/sei/base_sei.py | 34 +++++++++++++++++++ .../models/submodels/interface/sei/no_sei.py | 6 +++- .../interface/total_interfacial_current.py | 5 ++- src/pybamm/parameters/geometric_parameters.py | 6 ++-- src/pybamm/parameters/parameter_values.py | 1 - tests/unit/test_citations.py | 4 +-- .../test_base_battery_model.py | 4 +++ .../test_parameters/test_parameter_values.py | 2 +- 9 files changed, 66 insertions(+), 12 deletions(-) diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index c511d4ded0..71a68d5c46 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -241,6 +241,11 @@ def _get_standard_interfacial_current_variables(self, j): f"X-averaged {domain} electrode {reaction_name}" "interfacial current density [A.m-2]": j_av, } + if self.phase_name == reaction_name: + variables.update({ + f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]": j, + f"X-averaged {domain} electrode {reaction_name}interfacial current density [A.m-2]": j_av, + }) return variables @@ -309,14 +314,17 @@ def _get_standard_volumetric_current_density_variables(self, variables): a = variables[ f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]" ] - - if False: + if "plating" in reaction_name: + j = variables[ + f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]" + ] + elif phase_name in reaction_name: j = variables[ - f"{Domain} electrode {phase_name}{reaction_name}interfacial current density distribution [A.m-2]" + f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]" ] else: j = variables[ - f"{Domain} electrode {phase_name}interfacial current density [A.m-2]" + f"{Domain} electrode {phase_name}{reaction_name}interfacial current density [A.m-2]" ] print(reaction_name) a_j = a * j diff --git a/src/pybamm/models/submodels/interface/sei/base_sei.py b/src/pybamm/models/submodels/interface/sei/base_sei.py index ee8d6df7e3..56d1e228ca 100644 --- a/src/pybamm/models/submodels/interface/sei/base_sei.py +++ b/src/pybamm/models/submodels/interface/sei/base_sei.py @@ -263,3 +263,37 @@ def _get_standard_reaction_variables(self, j_sei): ) return variables + + def _get_standard_reaction_distribution_variables(self, j_sei): + """ + A private function to obtain the standard variables which + can be derived from the SEI interfacial reaction current + + Parameters + ---------- + j_sei : :class:`pybamm.Symbol` + The SEI interfacial reaction current. + + Returns + ------- + variables : dict + The variables which can be derived from the SEI currents. + """ + domain, Domain = self.domain_Domain + j_sei = pybamm.size_average(j_sei) + variables = { + f"{Domain} electrode {self.reaction_name}" + "interfacial current density [A.m-2]": j_sei, + } + + if self.reaction_loc != "interface": + j_sei_av = pybamm.x_average(j_sei) + variables.update( + { + f"X-averaged {domain} electrode {self.reaction_name}" + "interfacial current density [A.m-2]": j_sei_av, + } + ) + + return variables + diff --git a/src/pybamm/models/submodels/interface/sei/no_sei.py b/src/pybamm/models/submodels/interface/sei/no_sei.py index 242c8676ea..83eacb41e6 100644 --- a/src/pybamm/models/submodels/interface/sei/no_sei.py +++ b/src/pybamm/models/submodels/interface/sei/no_sei.py @@ -43,11 +43,15 @@ def get_fundamental_variables(self): pybamm.Scalar(0), f"{domain} electrode", "current collector" ) variables = self._get_standard_thickness_variables(zero) - variables.update(self._get_standard_reaction_variables(zero)) + if self.size_distribution: + variables.update(self._get_standard_reaction_distribution_variables(zero)) variables.update( self._get_standard_size_distribution_interfacial_current_variables(zero) ) + else: + variables.update(self._get_standard_reaction_variables(zero)) + variables.update(self._get_standard_interfacial_current_variables(zero)) return variables def get_coupled_variables(self, variables): diff --git a/src/pybamm/models/submodels/interface/total_interfacial_current.py b/src/pybamm/models/submodels/interface/total_interfacial_current.py index 79e13e37f6..36786307cf 100644 --- a/src/pybamm/models/submodels/interface/total_interfacial_current.py +++ b/src/pybamm/models/submodels/interface/total_interfacial_current.py @@ -165,7 +165,10 @@ def _get_whole_cell_coupled_variables(self, variables): if domain == "separator": var_dict[domain] = zero_s else: - var_dict[domain] = variables[variable_template.format(domain + " ")] + var = variables[variable_template.format(domain + " ")] + if "size" in var.domain[0]: + var = pybamm.size_average(var) + var_dict[domain] = var var = pybamm.concatenation(*var_dict.values()) variables.update({variable_template.format(""): var}) diff --git a/src/pybamm/parameters/geometric_parameters.py b/src/pybamm/parameters/geometric_parameters.py index 415ffb8e72..2998e97318 100644 --- a/src/pybamm/parameters/geometric_parameters.py +++ b/src/pybamm/parameters/geometric_parameters.py @@ -165,12 +165,14 @@ def f_a_dist(self, R): Dimensional electrode area-weighted particle-size distribution """ Domain = self.domain.capitalize() - if "primary" in R.domains["primary"][0]: + if len(R.domain) == 0: + self.phase_prefactor = "" + elif "primary" in R.domains["primary"][0]: self.phase_prefactor = "Primary: " elif "secondary" in R.domains["primary"][0]: self.phase_prefactor = "Secondary: " else: - pass + self.phase_prefactor = "" inputs = {f"{self.phase_prefactor}{Domain} particle-size variable [m]": R} return pybamm.FunctionParameter( f"{self.phase_prefactor}{Domain} " diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py index 318184c356..9ddbfd2fd3 100644 --- a/src/pybamm/parameters/parameter_values.py +++ b/src/pybamm/parameters/parameter_values.py @@ -666,7 +666,6 @@ def _process_symbol(self, symbol): elif isinstance(symbol, pybamm.FunctionParameter): function_name = self[symbol.name] - print(symbol.name) if isinstance( function_name, (numbers.Number, pybamm.Interpolant, pybamm.InputParameter), diff --git a/tests/unit/test_citations.py b/tests/unit/test_citations.py index 9b1c9fa2d3..dccef9b51c 100644 --- a/tests/unit/test_citations.py +++ b/tests/unit/test_citations.py @@ -357,13 +357,13 @@ def test_sripad_2020(self): citations._reset() assert "Sripad2020" not in citations._papers_to_cite - pybamm.kinetics.Marcus(None, None, None, None, None) + pybamm.kinetics.Marcus(None, "negative", None, None, None) assert "Sripad2020" in citations._papers_to_cite assert "Sripad2020" in citations._citation_tags.keys() citations._reset() assert "Sripad2020" not in citations._papers_to_cite - pybamm.kinetics.MarcusHushChidsey(None, None, None, None, None) + pybamm.kinetics.MarcusHushChidsey(None, "negative", None, None, None) assert "Sripad2020" in citations._papers_to_cite assert "Sripad2020" in citations._citation_tags.keys() diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py index be7ee861c4..554b65d880 100644 --- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py +++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py @@ -140,6 +140,10 @@ def test_default_var_pts(self): "z": 10, "R_n": 30, "R_p": 30, + "R_n_prim": 30, + "R_p_prim": 30, + "R_n_sec": 30, + "R_p_sec": 30, } model = pybamm.BaseBatteryModel({"dimensionality": 0}) assert var_pts == model.default_var_pts diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index cdcfb30ede..0e32f60476 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -893,7 +893,7 @@ def dist(R): var_av_proc = param.process_symbol(var_av) assert isinstance(var_av_proc, pybamm.SizeAverage) - R = pybamm.SpatialVariable("R", "negative particle size") + R = pybamm.SpatialVariable("R_n", "negative particle size") assert var_av_proc.f_a_dist == R**2 def test_process_not_constant(self): From 19c1febb87e537a7a64a31247bb6997bd4b83ed7 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Tue, 17 Dec 2024 19:08:58 -0800 Subject: [PATCH 07/17] add test --- .../submodels/interface/base_interface.py | 15 +++++--- .../base_lithium_ion_tests.py | 35 +++++++++++++++++++ .../test_lithium_ion/test_newman_tobias.py | 3 ++ 3 files changed, 48 insertions(+), 5 deletions(-) diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index 71a68d5c46..f14216787e 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -232,7 +232,10 @@ def _get_standard_interfacial_current_variables(self, j): j_av = pybamm.x_average(j) if j.domain == []: j = pybamm.FullBroadcast(j, f"{domain} electrode", "current collector") - elif j.domain == ["current collector"]: + elif ( + j.domain == ["current collector"] + and getattr(self, "reaction_loc", None) != "interface" + ): j = pybamm.PrimaryBroadcast(j, f"{domain} electrode") variables = { @@ -242,10 +245,12 @@ def _get_standard_interfacial_current_variables(self, j): "interfacial current density [A.m-2]": j_av, } if self.phase_name == reaction_name: - variables.update({ - f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]": j, - f"X-averaged {domain} electrode {reaction_name}interfacial current density [A.m-2]": j_av, - }) + variables.update( + { + f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]": j, + f"X-averaged {domain} electrode {reaction_name}interfacial current density [A.m-2]": j_av, + } + ) return variables diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py index 768872e454..5f0d4d7d83 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py @@ -351,6 +351,41 @@ def test_composite_graphite_silicon(self): ) self.run_basic_processing_test(options, parameter_values=parameter_values) + def test_particle_size_composite(self): + options = { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + "particle size": "distribution", + "surface form": "algebraic", + } + parameter_values = pybamm.ParameterValues("Chen2020_composite") + name = "Negative electrode active material volume fraction" + x = 0.1 + parameter_values.update( + {f"Primary: {name}": (1 - x) * 0.75, f"Secondary: {name}": x * 0.75} + ) + parameter_values = pybamm.get_size_distribution_parameters( + parameter_values, + composite="negative", + R_min_n_prim=0.9, + R_min_n_sec=0.9, + R_max_n_prim=1.1, + R_max_n_sec=1.1, + ) + # self.run_basic_processing_test(options, parameter_values=parameter_values) + # TODO: fix this test + sim = pybamm.Simulation(self.model(options), parameter_values=parameter_values) + sim.solve(t_eval=np.linspace(0, 3600)) + + def test_particle_size_distribution(self): + options = { + "particle size": "distribution", + "surface form": "algebraic", + } + parameter_values = self.model(options).default_parameter_values + parameter_values = pybamm.get_size_distribution_parameters(parameter_values) + self.run_basic_processing_test(options, parameter_values=parameter_values) + def test_composite_graphite_silicon_sei(self): options = { "particle phases": ("2", "1"), diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py index fb091ef56d..c24a063f70 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py @@ -31,6 +31,9 @@ def test_composite_graphite_silicon_sei(self): def test_composite_reaction_driven_LAM(self): pass # skip this test + def test_particle_size_composite(self): + pass # skip this test + def test_composite_stress_driven_LAM(self): pass # skip this test From 03a367cd8eab440de6ad1a4d7ca81b6467983dad Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 18 Dec 2024 03:10:41 +0000 Subject: [PATCH 08/17] style: pre-commit fixes --- src/pybamm/models/submodels/interface/sei/base_sei.py | 3 +-- src/pybamm/models/submodels/interface/sei/no_sei.py | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/pybamm/models/submodels/interface/sei/base_sei.py b/src/pybamm/models/submodels/interface/sei/base_sei.py index 56d1e228ca..d4a3711c1b 100644 --- a/src/pybamm/models/submodels/interface/sei/base_sei.py +++ b/src/pybamm/models/submodels/interface/sei/base_sei.py @@ -263,7 +263,7 @@ def _get_standard_reaction_variables(self, j_sei): ) return variables - + def _get_standard_reaction_distribution_variables(self, j_sei): """ A private function to obtain the standard variables which @@ -296,4 +296,3 @@ def _get_standard_reaction_distribution_variables(self, j_sei): ) return variables - diff --git a/src/pybamm/models/submodels/interface/sei/no_sei.py b/src/pybamm/models/submodels/interface/sei/no_sei.py index 83eacb41e6..a0cfec6ccf 100644 --- a/src/pybamm/models/submodels/interface/sei/no_sei.py +++ b/src/pybamm/models/submodels/interface/sei/no_sei.py @@ -43,7 +43,7 @@ def get_fundamental_variables(self): pybamm.Scalar(0), f"{domain} electrode", "current collector" ) variables = self._get_standard_thickness_variables(zero) - + if self.size_distribution: variables.update(self._get_standard_reaction_distribution_variables(zero)) variables.update( From 55dd848e246b2638e536819dbdecfd3ab455354a Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Tue, 17 Dec 2024 19:16:06 -0800 Subject: [PATCH 09/17] cleanup and changelog --- CHANGELOG.md | 1 + src/pybamm/expression_tree/averages.py | 1 - src/pybamm/meshes/meshes.py | 1 - src/pybamm/models/submodels/interface/base_interface.py | 1 - 4 files changed, 1 insertion(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ebc582547b..a5179ae4ee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- Made composite electrode model compatible with particle size distribution ([#4687](https://github.com/pybamm-team/PyBaMM/pull/4687)) - Added `Symbol.post_order()` method to return an iterable that steps through the tree in post-order fashion. ([#4684](https://github.com/pybamm-team/PyBaMM/pull/4684)) - Added two more submodels (options) for the SEI: Lars von Kolzenberg (2020) model and Tunneling Limit model ([#4394](https://github.com/pybamm-team/PyBaMM/pull/4394)) diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py index 33f38c68b3..a14683c873 100644 --- a/src/pybamm/expression_tree/averages.py +++ b/src/pybamm/expression_tree/averages.py @@ -388,7 +388,6 @@ def size_average( ] for domain in list(symbol.domains.values()) ): - print(f"symbol: {symbol.domain}") return symbol # If symbol is a primary broadcast to "particle size", take the orphan diff --git a/src/pybamm/meshes/meshes.py b/src/pybamm/meshes/meshes.py index 45006bf88a..3ec291b1b8 100644 --- a/src/pybamm/meshes/meshes.py +++ b/src/pybamm/meshes/meshes.py @@ -178,7 +178,6 @@ def combine_submeshes(self, *submeshnames): raise pybamm.DomainError( "trying to combine two meshes in different coordinate systems" ) - print(submeshnames) combined_submesh_edges = np.concatenate( [self[submeshnames[0]].edges] + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]] diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index f14216787e..337ac1ce86 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -331,7 +331,6 @@ def _get_standard_volumetric_current_density_variables(self, variables): j = variables[ f"{Domain} electrode {phase_name}{reaction_name}interfacial current density [A.m-2]" ] - print(reaction_name) a_j = a * j a_j_av = pybamm.x_average(a_j) From 29f4f3ddfd5a8d0237e6dcea882f9c80355a98f7 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Wed, 18 Dec 2024 08:06:16 -0800 Subject: [PATCH 10/17] fix notebook issue --- docs/source/examples/notebooks/models/using-submodels.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/examples/notebooks/models/using-submodels.ipynb b/docs/source/examples/notebooks/models/using-submodels.ipynb index cce32fcefd..db889972f8 100644 --- a/docs/source/examples/notebooks/models/using-submodels.ipynb +++ b/docs/source/examples/notebooks/models/using-submodels.ipynb @@ -527,10 +527,10 @@ " model.param, \"positive\", model.options, cracks=True\n", ")\n", "model.submodels[\"Negative lithium plating\"] = pybamm.lithium_plating.NoPlating(\n", - " model.param, \"Negative\"\n", + " model.param, \"negative\"\n", ")\n", "model.submodels[\"Positive lithium plating\"] = pybamm.lithium_plating.NoPlating(\n", - " model.param, \"Positive\"\n", + " model.param, \"positive\"\n", ")" ] }, From 835eac65653bf300901123bdde2871dce610f121 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Wed, 18 Dec 2024 15:23:33 -0800 Subject: [PATCH 11/17] valentin style changes --- src/pybamm/expression_tree/averages.py | 27 +++++++------------ .../submodels/interface/sei/base_sei.py | 3 +-- 2 files changed, 10 insertions(+), 20 deletions(-) diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py index a14683c873..841cb37751 100644 --- a/src/pybamm/expression_tree/averages.py +++ b/src/pybamm/expression_tree/averages.py @@ -408,30 +408,21 @@ def size_average( geo = pybamm.geometric_parameters if "negative" in symbol.domain[0]: if "primary" in symbol.domain[0]: - R = pybamm.SpatialVariable( - "R_n_prim", domains=symbol.domains, coord_sys="cartesian" - ) + name = "R_n_prim" elif "secondary" in symbol.domain[0]: - R = pybamm.SpatialVariable( - "R_n_sec", domains=symbol.domains, coord_sys="cartesian" - ) + name = "R_n_sec" else: - R = pybamm.SpatialVariable( - "R_n", domains=symbol.domains, coord_sys="cartesian" - ) + name = "R_n" else: if "primary" in symbol.domain[0]: - R = pybamm.SpatialVariable( - "R_p_prim", domains=symbol.domains, coord_sys="cartesian" - ) + name = "R_p_prim" elif "secondary" in symbol.domain[0]: - R = pybamm.SpatialVariable( - "R_p_sec", domains=symbol.domains, coord_sys="cartesian" - ) + name = "R_p_sec" else: - R = pybamm.SpatialVariable( - "R_p", domains=symbol.domains, coord_sys="cartesian" - ) + name = "R_p" + R = pybamm.SpatialVariable( + name, domains=symbol.domains, coord_sys="cartesian" + ) if ["negative particle size"] in symbol.domains.values() or [ "negative primary particle size" ] in symbol.domains.values(): diff --git a/src/pybamm/models/submodels/interface/sei/base_sei.py b/src/pybamm/models/submodels/interface/sei/base_sei.py index d4a3711c1b..af5a929bb5 100644 --- a/src/pybamm/models/submodels/interface/sei/base_sei.py +++ b/src/pybamm/models/submodels/interface/sei/base_sei.py @@ -155,10 +155,9 @@ def _get_standard_concentration_variables(self, variables): if self.size_distribution: n_SEI_sav = pybamm.size_average(n_SEI) n_SEI_xav = pybamm.x_average(n_SEI_sav) - n_SEI_av = pybamm.yz_average(n_SEI_xav) else: n_SEI_xav = pybamm.x_average(n_SEI) - n_SEI_av = pybamm.yz_average(n_SEI_xav) + n_SEI_av = pybamm.yz_average(n_SEI_xav) # Calculate change in SEI concentration with respect to initial state delta_n_SEI = n_SEI_av - n_SEI_0 From fc99bbb289e7b05a17e709f7b0280c5c09397504 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Wed, 18 Dec 2024 17:23:31 -0800 Subject: [PATCH 12/17] fully test distribution-parameters --- .../submodels/interface/base_interface.py | 9 +- .../size_distribution_parameters.py | 138 ++++++++++++------ .../test_models/standard_output_tests.py | 8 +- .../base_lithium_ion_tests.py | 35 ----- .../test_lithium_ion/test_dfn.py | 29 ++++ .../test_size_distribution_parameters.py | 59 ++++++++ 6 files changed, 193 insertions(+), 85 deletions(-) diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index 337ac1ce86..7d356313ab 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -364,7 +364,14 @@ def _get_standard_overpotential_variables(self, eta_r): # Size average. For eta_r variables that depend on particle size, see # "_get_standard_size_distribution_overpotential_variables" - if eta_r.domain in [["negative particle size"], ["positive particle size"]]: + if eta_r.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: eta_r = pybamm.size_average(eta_r) # X-average, and broadcast if necessary diff --git a/src/pybamm/parameters/size_distribution_parameters.py b/src/pybamm/parameters/size_distribution_parameters.py index 6693cd0af3..6b7087b8c8 100644 --- a/src/pybamm/parameters/size_distribution_parameters.py +++ b/src/pybamm/parameters/size_distribution_parameters.py @@ -12,18 +12,26 @@ def get_size_distribution_parameters( R_n_av_prim=None, R_n_av_sec=None, R_p_av=None, + R_p_av_prim=None, + R_p_av_sec=None, sd_n=0.3, sd_n_prim=0.3, sd_n_sec=0.3, sd_p=0.3, + sd_p_prim=0.3, + sd_p_sec=0.3, R_min_n=None, R_min_n_prim=None, R_min_n_sec=None, R_min_p=None, + R_min_p_prim=None, + R_min_p_sec=None, R_max_n=None, R_max_n_prim=None, R_max_n_sec=None, R_max_p=None, + R_max_p_prim=None, + R_max_p_sec=None, working_electrode="both", composite="neither", ): @@ -71,30 +79,7 @@ def get_size_distribution_parameters( """ if working_electrode == "both": # Radii from given parameter set - if composite != "negative" and composite != "both": - R_n_typ = param["Negative particle radius [m]"] - # Set the mean particle radii - R_n_av = R_n_av or R_n_typ - # Minimum radii - R_min_n = R_min_n or np.max([0, 1 - sd_n * 5]) - - # Max radii - R_max_n = R_max_n or (1 + sd_n * 5) - - # Area-weighted particle-size distribution - def f_a_dist_n(R): - return lognormal(R, R_n_av, sd_n * R_n_av) - - param.update( - { - "Negative minimum particle radius [m]": R_min_n * R_n_av, - "Negative maximum particle radius [m]": R_max_n * R_n_av, - "Negative area-weighted " - + "particle-size distribution [m-1]": f_a_dist_n, - }, - check_already_exists=False, - ) - else: + if composite == "negative" or composite == "both": R_n_typ_prim = param["Primary: Negative particle radius [m]"] R_n_typ_sec = param["Secondary: Negative particle radius [m]"] # Set the mean particle radii @@ -132,31 +117,94 @@ def f_a_dist_n_sec(R): }, check_already_exists=False, ) + else: + R_n_typ = param["Negative particle radius [m]"] + # Set the mean particle radii + R_n_av = R_n_av or R_n_typ + # Minimum radii + R_min_n = R_min_n or np.max([0, 1 - sd_n * 5]) - # Radii from given parameter set - R_p_typ = param["Positive particle radius [m]"] - - # Set the mean particle radii - R_p_av = R_p_av or R_p_typ - - # Minimum radii - R_min_p = R_min_p or np.max([0, 1 - sd_p * 5]) + # Max radii + R_max_n = R_max_n or (1 + sd_n * 5) - # Max radii - R_max_p = R_max_p or (1 + sd_p * 5) + # Area-weighted particle-size distribution + def f_a_dist_n(R): + return lognormal(R, R_n_av, sd_n * R_n_av) - # Area-weighted particle-size distribution - def f_a_dist_p(R): - return lognormal(R, R_p_av, sd_p * R_p_av) + param.update( + { + "Negative minimum particle radius [m]": R_min_n * R_n_av, + "Negative maximum particle radius [m]": R_max_n * R_n_av, + "Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n, + }, + check_already_exists=False, + ) - param.update( - { - "Positive minimum particle radius [m]": R_min_p * R_p_av, - "Positive maximum particle radius [m]": R_max_p * R_p_av, - "Positive area-weighted " + "particle-size distribution [m-1]": f_a_dist_p, - }, - check_already_exists=False, - ) + if composite == "positive" or composite == "both": + R_p_typ_prim = param["Primary: Positive particle radius [m]"] + R_p_typ_sec = param["Secondary: Positive particle radius [m]"] + # Set the mean particle radii + R_p_av_prim = R_p_av_prim or R_p_typ_prim + R_p_av_sec = R_p_av_sec or R_p_typ_sec + # Minimum radii + R_min_p_prim = R_min_p_prim or np.max([0, 1 - sd_p_prim * 5]) + R_min_p_sec = R_min_p_sec or np.max([0, 1 - sd_p_sec * 5]) + + # Max radii + R_max_p_prim = R_max_p_prim or (1 + sd_p_prim * 5) + R_max_p_sec = R_max_p_sec or (1 + sd_p_sec * 5) + + # Area-weighted particle-size distribution + def f_a_dist_p_prim(R): + return lognormal(R, R_p_av_prim, sd_p_prim * R_p_av_prim) + + def f_a_dist_p_sec(R): + return lognormal(R, R_p_av_sec, sd_p_sec * R_p_av_sec) + + param.update( + { + "Primary: Positive minimum particle radius [m]": R_min_p_prim + * R_p_av_prim, + "Primary: Positive maximum particle radius [m]": R_max_p_prim + * R_p_av_prim, + "Primary: Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p_prim, + "Secondary: Positive minimum particle radius [m]": R_min_p_sec + * R_p_av_sec, + "Secondary: Positive maximum particle radius [m]": R_max_p_sec + * R_p_av_sec, + "Secondary: Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p_sec, + }, + check_already_exists=False, + ) + else: + # Radii from given parameter set + R_p_typ = param["Positive particle radius [m]"] + + # Set the mean particle radii + R_p_av = R_p_av or R_p_typ + + # Minimum radii + R_min_p = R_min_p or np.max([0, 1 - sd_p * 5]) + + # Max radii + R_max_p = R_max_p or (1 + sd_p * 5) + + # Area-weighted particle-size distribution + def f_a_dist_p(R): + return lognormal(R, R_p_av, sd_p * R_p_av) + + param.update( + { + "Positive minimum particle radius [m]": R_min_p * R_p_av, + "Positive maximum particle radius [m]": R_max_p * R_p_av, + "Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p, + }, + check_already_exists=False, + ) return param diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py index ca5e27607d..4e6788724f 100644 --- a/tests/integration/test_models/standard_output_tests.py +++ b/tests/integration/test_models/standard_output_tests.py @@ -322,17 +322,17 @@ def __init__(self, model, param, disc, solution, operating_condition): # Take only the x-averaged of these for now, since variables cannot have # 4 domains yet self.c_s_n_dist = solution[ - "X-averaged negative particle concentration distribution" + f"X-averaged negative {self.phase_name_n}particle concentration distribution" ] self.c_s_p_dist = solution[ - "X-averaged positive particle concentration distribution" + f"X-averaged positive {self.phase_name_p}particle concentration distribution" ] self.c_s_n_surf_dist = solution[ - "Negative particle surface concentration distribution" + f"Negative {self.phase_name_n}particle surface concentration distribution" ] self.c_s_p_surf_dist = solution[ - "Positive particle surface concentration distribution" + f"Positive {self.phase_name_p}particle surface concentration distribution" ] def test_concentration_increase_decrease(self): diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py index 5f0d4d7d83..768872e454 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py @@ -351,41 +351,6 @@ def test_composite_graphite_silicon(self): ) self.run_basic_processing_test(options, parameter_values=parameter_values) - def test_particle_size_composite(self): - options = { - "particle phases": ("2", "1"), - "open-circuit potential": (("single", "current sigmoid"), "single"), - "particle size": "distribution", - "surface form": "algebraic", - } - parameter_values = pybamm.ParameterValues("Chen2020_composite") - name = "Negative electrode active material volume fraction" - x = 0.1 - parameter_values.update( - {f"Primary: {name}": (1 - x) * 0.75, f"Secondary: {name}": x * 0.75} - ) - parameter_values = pybamm.get_size_distribution_parameters( - parameter_values, - composite="negative", - R_min_n_prim=0.9, - R_min_n_sec=0.9, - R_max_n_prim=1.1, - R_max_n_sec=1.1, - ) - # self.run_basic_processing_test(options, parameter_values=parameter_values) - # TODO: fix this test - sim = pybamm.Simulation(self.model(options), parameter_values=parameter_values) - sim.solve(t_eval=np.linspace(0, 3600)) - - def test_particle_size_distribution(self): - options = { - "particle size": "distribution", - "surface form": "algebraic", - } - parameter_values = self.model(options).default_parameter_values - parameter_values = pybamm.get_size_distribution_parameters(parameter_values) - self.run_basic_processing_test(options, parameter_values=parameter_values) - def test_composite_graphite_silicon_sei(self): options = { "particle phases": ("2", "1"), diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py index 5ede55ff6d..7a728a744b 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py @@ -52,6 +52,10 @@ def setup(self): "R_p": 3, "y": 5, "z": 5, + "R_n_prim": 3, + "R_n_sec": 3, + "R_p_prim": 3, + "R_p_sec": 3, } def test_basic_processing(self): @@ -62,6 +66,31 @@ def test_basic_processing(self): ) modeltest.test_all() + def test_composite(self): + options = { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + "particle size": "distribution", + } + parameter_values = pybamm.ParameterValues("Chen2020_composite") + name = "Negative electrode active material volume fraction" + x = 0.1 + parameter_values.update( + {f"Primary: {name}": (1 - x) * 0.75, f"Secondary: {name}": x * 0.75} + ) + parameter_values = pybamm.get_size_distribution_parameters( + parameter_values, + composite="negative", + R_min_n_prim=0.9, + R_min_n_sec=0.9, + R_max_n_prim=1.1, + R_max_n_sec=1.1, + ) + # self.run_basic_processing_test(options, parameter_values=parameter_values) + model = pybamm.lithium_ion.DFN(options) + modeltest = tests.StandardModelTest(model, parameter_values=parameter_values) + modeltest.test_all() + def test_basic_processing_tuple(self): options = {"particle size": ("single", "distribution")} model = pybamm.lithium_ion.DFN(options) diff --git a/tests/unit/test_parameters/test_size_distribution_parameters.py b/tests/unit/test_parameters/test_size_distribution_parameters.py index 414b422055..a90c6a3390 100644 --- a/tests/unit/test_parameters/test_size_distribution_parameters.py +++ b/tests/unit/test_parameters/test_size_distribution_parameters.py @@ -40,3 +40,62 @@ def test_parameter_values(self): R_test = pybamm.Scalar(1.0) values.evaluate(param.n.prim.f_a_dist(R_test)) values.evaluate(param.p.prim.f_a_dist(R_test)) + + # check that the composite parameters works as well + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="negative", + ) + assert "Primary: Negative maximum particle radius [m]" in values_composite + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="positive", + ) + assert "Primary: Positive maximum particle radius [m]" in values_composite + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="both", + ) + assert "Primary: Negative maximum particle radius [m]" in values_composite + assert "Primary: Positive maximum particle radius [m]" in values_composite + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + ) + assert "Primary: Negative maximum particle radius [m]" not in values_composite + assert "Primary: Positive maximum particle radius [m]" not in values_composite From cd9fd3603545049b129e9cc9dfdb5ff19c03372b Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Wed, 18 Dec 2024 17:59:18 -0800 Subject: [PATCH 13/17] coverage --- tests/unit/test_expression_tree/test_averages.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/tests/unit/test_expression_tree/test_averages.py b/tests/unit/test_expression_tree/test_averages.py index 5f9736693e..9065b85ca9 100644 --- a/tests/unit/test_expression_tree/test_averages.py +++ b/tests/unit/test_expression_tree/test_averages.py @@ -201,7 +201,14 @@ def test_size_average(self): ) assert average_a == a - for domain in [["negative particle size"], ["positive particle size"]]: + for domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: a = pybamm.Symbol("a", domain=domain) R = pybamm.SpatialVariable("R", domain) av_a = pybamm.size_average(a) From 505495eeb6032b8f21500acb75049d1cd5e50d77 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Wed, 18 Dec 2024 18:02:55 -0800 Subject: [PATCH 14/17] more coverage --- .../test_full_battery_models/test_lithium_ion/test_dfn.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py index 973a0f348b..56f4c76f44 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py @@ -27,6 +27,10 @@ def test_well_posed_size_distribution_tuple(self): options = {"particle size": ("single", "distribution")} self.check_well_posedness(options) + def test_well_posed_size_distribution_composite(self): + options = {"particle size": "distribution", "particle phases": "2"} + self.check_well_posedness(options) + def test_well_posed_current_sigmoid_ocp_with_psd(self): options = { "open-circuit potential": "current sigmoid", From b33107480547a879cd4f060047b84a0969229eaa Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Wed, 18 Dec 2024 18:47:04 -0800 Subject: [PATCH 15/17] remove untouched lines --- .../models/submodels/interface/base_interface.py | 15 +++------------ .../interface/total_interfacial_current.py | 2 -- tests/unit/test_geometry/test_battery_geometry.py | 10 ++++++++++ .../test_size_distribution_parameters.py | 4 ++++ 4 files changed, 17 insertions(+), 14 deletions(-) diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index 7d356313ab..83f7c81b43 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -319,18 +319,9 @@ def _get_standard_volumetric_current_density_variables(self, variables): a = variables[ f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]" ] - if "plating" in reaction_name: - j = variables[ - f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]" - ] - elif phase_name in reaction_name: - j = variables[ - f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]" - ] - else: - j = variables[ - f"{Domain} electrode {phase_name}{reaction_name}interfacial current density [A.m-2]" - ] + j = variables[ + f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]" + ] a_j = a * j a_j_av = pybamm.x_average(a_j) diff --git a/src/pybamm/models/submodels/interface/total_interfacial_current.py b/src/pybamm/models/submodels/interface/total_interfacial_current.py index 36786307cf..3df4ff757a 100644 --- a/src/pybamm/models/submodels/interface/total_interfacial_current.py +++ b/src/pybamm/models/submodels/interface/total_interfacial_current.py @@ -166,8 +166,6 @@ def _get_whole_cell_coupled_variables(self, variables): var_dict[domain] = zero_s else: var = variables[variable_template.format(domain + " ")] - if "size" in var.domain[0]: - var = pybamm.size_average(var) var_dict[domain] = var var = pybamm.concatenation(*var_dict.values()) variables.update({variable_template.format(""): var}) diff --git a/tests/unit/test_geometry/test_battery_geometry.py b/tests/unit/test_geometry/test_battery_geometry.py index c73ef89af1..30d417fcbd 100644 --- a/tests/unit/test_geometry/test_battery_geometry.py +++ b/tests/unit/test_geometry/test_battery_geometry.py @@ -77,6 +77,16 @@ def test_geometry(self): geometry["positive secondary particle"]["r_p_sec"]["max"] == geo.p.sec.R_typ ) + options = { + "particle phases": "2", + "particle size": "distribution", + } + geometry = pybamm.battery_geometry(options=options) + assert "negative primary particle size" in geometry + assert "positive primary particle size" in geometry + assert "negative secondary particle size" in geometry + assert "positive secondary particle size" in geometry + def test_geometry_error(self): with pytest.raises(pybamm.GeometryError, match="Invalid current"): pybamm.battery_geometry( diff --git a/tests/unit/test_parameters/test_size_distribution_parameters.py b/tests/unit/test_parameters/test_size_distribution_parameters.py index a90c6a3390..5b6cef0aaa 100644 --- a/tests/unit/test_parameters/test_size_distribution_parameters.py +++ b/tests/unit/test_parameters/test_size_distribution_parameters.py @@ -56,6 +56,8 @@ def test_parameter_values(self): composite="negative", ) assert "Primary: Negative maximum particle radius [m]" in values_composite + values.evaluate(param.n.prim.f_a_dist_n_prim(R_test)) + values.evaluate(param.n.sec.f_a_dist_n_sec(R_test)) params_composite = pybamm.ParameterValues("Chen2020_composite") params_composite.update( { @@ -70,6 +72,8 @@ def test_parameter_values(self): composite="positive", ) assert "Primary: Positive maximum particle radius [m]" in values_composite + values.evaluate(param.p.prim.f_a_dist_p_prim(R_test)) + values.evaluate(param.p.sec.f_a_dist_p_sec(R_test)) params_composite = pybamm.ParameterValues("Chen2020_composite") params_composite.update( { From 84fd34207f30f72127e2f56dd6a29f6fca5a7892 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Wed, 18 Dec 2024 19:16:17 -0800 Subject: [PATCH 16/17] fix failing test --- .../test_size_distribution_parameters.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/tests/unit/test_parameters/test_size_distribution_parameters.py b/tests/unit/test_parameters/test_size_distribution_parameters.py index 5b6cef0aaa..a0f7847a32 100644 --- a/tests/unit/test_parameters/test_size_distribution_parameters.py +++ b/tests/unit/test_parameters/test_size_distribution_parameters.py @@ -56,8 +56,13 @@ def test_parameter_values(self): composite="negative", ) assert "Primary: Negative maximum particle radius [m]" in values_composite - values.evaluate(param.n.prim.f_a_dist_n_prim(R_test)) - values.evaluate(param.n.sec.f_a_dist_n_sec(R_test)) + param = pybamm.LithiumIonParameters({"particle phases": ("2", "1")}) + R_test_n_prim = pybamm.Scalar(1e-6) + R_test_n_prim.domains = {"primary": ["negative primary particle size"]} + R_test_n_sec = pybamm.Scalar(1e-6) + R_test_n_sec.domains = {"primary": ["negative secondary particle size"]} + values_composite.evaluate(param.n.prim.f_a_dist(R_test_n_prim)) + values_composite.evaluate(param.n.sec.f_a_dist(R_test_n_sec)) params_composite = pybamm.ParameterValues("Chen2020_composite") params_composite.update( { @@ -72,8 +77,13 @@ def test_parameter_values(self): composite="positive", ) assert "Primary: Positive maximum particle radius [m]" in values_composite - values.evaluate(param.p.prim.f_a_dist_p_prim(R_test)) - values.evaluate(param.p.sec.f_a_dist_p_sec(R_test)) + param = pybamm.LithiumIonParameters({"particle phases": ("1", "2")}) + R_test_p_prim = pybamm.Scalar(1e-6) + R_test_p_prim.domains = {"primary": ["positive primary particle size"]} + R_test_p_sec = pybamm.Scalar(1e-6) + R_test_p_sec.domains = {"primary": ["positive secondary particle size"]} + values_composite.evaluate(param.p.prim.f_a_dist(R_test_p_prim)) + values_composite.evaluate(param.p.sec.f_a_dist(R_test_p_sec)) params_composite = pybamm.ParameterValues("Chen2020_composite") params_composite.update( { From f545b5ce9d07bb9a9c463ef1d6ae38baa57aa13a Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Thu, 19 Dec 2024 09:16:16 -0800 Subject: [PATCH 17/17] rob comments --- src/pybamm/expression_tree/averages.py | 21 +++++++------------ .../interface/lithium_plating/base_plating.py | 4 ---- .../interface/lithium_plating/no_plating.py | 8 +++---- .../models/submodels/interface/sei/no_sei.py | 2 +- 4 files changed, 13 insertions(+), 22 deletions(-) diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py index 841cb37751..c95ccb4e5e 100644 --- a/src/pybamm/expression_tree/averages.py +++ b/src/pybamm/expression_tree/averages.py @@ -406,20 +406,15 @@ def size_average( else: if f_a_dist is None: geo = pybamm.geometric_parameters + name = "R" if "negative" in symbol.domain[0]: - if "primary" in symbol.domain[0]: - name = "R_n_prim" - elif "secondary" in symbol.domain[0]: - name = "R_n_sec" - else: - name = "R_n" - else: - if "primary" in symbol.domain[0]: - name = "R_p_prim" - elif "secondary" in symbol.domain[0]: - name = "R_p_sec" - else: - name = "R_p" + name += "_n" + elif "positive" in symbol.domain[0]: + name += "_p" + if "primary" in symbol.domain[0]: + name += "_prim" + elif "secondary" in symbol.domain[0]: + name += "_sec" R = pybamm.SpatialVariable( name, domains=symbol.domains, coord_sys="cartesian" ) diff --git a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py index 9425e375b2..c34eca9e05 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py @@ -90,10 +90,6 @@ def _get_standard_concentration_variables(self, c_plated_Li, c_dead_Li): L_dead_Li_av = pybamm.x_average(L_dead_Li) Q_dead_Li = c_dead_Li_av * L_k * self.param.L_y * self.param.L_z - # if self.size_distribution is True: - # Q_plated_Li = pybamm.size_average(Q_plated_Li) - # Q_dead_Li = pybamm.size_average(Q_dead_Li) - variables = { f"{Domain} {phase_name}lithium plating concentration " "[mol.m-3]": c_plated_Li, diff --git a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py index 64af5876df..8efa664a92 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py @@ -35,16 +35,16 @@ def get_fundamental_variables(self): }, ) variables = self._get_standard_concentration_variables(zero, zero) - if self.size_distribution is False: - variables.update(self._get_standard_overpotential_variables(zero)) - variables.update(self._get_standard_reaction_variables(zero)) - else: + if self.size_distribution: variables.update( self._get_standard_size_distribution_overpotential_variables(zero) ) variables.update( self._get_standard_size_distribution_reaction_variables(zero) ) + else: + variables.update(self._get_standard_reaction_variables(zero)) + variables.update(self._get_standard_overpotential_variables(zero)) return variables diff --git a/src/pybamm/models/submodels/interface/sei/no_sei.py b/src/pybamm/models/submodels/interface/sei/no_sei.py index a0cfec6ccf..10192a485f 100644 --- a/src/pybamm/models/submodels/interface/sei/no_sei.py +++ b/src/pybamm/models/submodels/interface/sei/no_sei.py @@ -51,7 +51,7 @@ def get_fundamental_variables(self): ) else: variables.update(self._get_standard_reaction_variables(zero)) - variables.update(self._get_standard_interfacial_current_variables(zero)) + variables.update(self._get_standard_interfacial_current_variables(zero)) return variables def get_coupled_variables(self, variables):