particle tuple pybamm-team#2669

pybamm/geometry/battery_geometry.py

@@ -63,14 +63,19 @@ def battery_geometry(
                     }
                 )
     # Add particle size domains
-    if options is not None and options["particle size"] == "distribution":
+    if options is not None and options.negative["particle size"] == "distribution":
         R_min_n = geo.n.prim.R_min
-        R_min_p = geo.p.prim.R_min
         R_max_n = geo.n.prim.R_max
-        R_max_p = geo.p.prim.R_max
         geometry.update(
             {
                 "negative particle size": {"R_n": {"min": R_min_n, "max": R_max_n}},
+            }
+        )
+    if options is not None and options.positive["particle size"] == "distribution":
+        R_min_p = geo.p.prim.R_min
+        R_max_p = geo.p.prim.R_max
+        geometry.update(
+            {
                 "positive particle size": {"R_p": {"min": R_min_p, "max": R_max_p}},
             }
         )

pybamm/models/full_battery_models/base_battery_model.py

@@ -73,11 +73,6 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "stress-driven", "reaction-driven", or "stress and reaction-driven".
                 A 2-tuple can be provided for different behaviour in negative and
                 positive electrodes.
-            * "particle phases": str
-                Number of phases present in the electrode. A 2-tuple can be provided for
-                different behaviour in negative and positive electrodes.
-                For example, set to ("2", "1") for a negative electrode with 2 phases,
-                e.g. graphite and silicon.
             * "operating mode" : str
                 Sets the operating mode for the model. This determines how the current
                 is set. Can be:
@@ -97,7 +92,19 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "particle" : str
                 Sets the submodel to use to describe behaviour within the particle.
                 Can be "Fickian diffusion" (default), "uniform profile",
-                "quadratic profile", or "quartic profile".
+                "quadratic profile", or "quartic profile". A 2-tuple can be provided for 
+                different behaviour in negative and positive electrodes.
+            * "particle mechanics" : str
+                Sets the model to account for mechanical effects such as particle
+                swelling and cracking. Can be "none" (default), "swelling only",
+                or "swelling and cracking".
+                A 2-tuple can be provided for different behaviour in negative and
+                positive electrodes.
+            * "particle phases": str
+                Number of phases present in the electrode. A 2-tuple can be provided for
+                different behaviour in negative and positive electrodes.
+                For example, set to ("2", "1") for a negative electrode with 2 phases,
+                e.g. graphite and silicon.
             * "particle shape" : str
                 Sets the model shape of the electrode particles. This is used to
                 calculate the surface area to volume ratio. Can be "spherical"
@@ -106,12 +113,6 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 Sets the model to include a single active particle size or a
                 distribution of sizes at any macroscale location. Can be "single"
                 (default) or "distribution". Option applies to both electrodes.
-            * "particle mechanics" : str
-                Sets the model to account for mechanical effects such as particle
-                swelling and cracking. Can be "none" (default), "swelling only",
-                or "swelling and cracking".
-                A 2-tuple can be provided for different behaviour in negative and
-                positive electrodes.
             * "SEI" : str
                 Set the SEI submodel to be used. Options are:
 
@@ -528,7 +529,7 @@ def __init__(self, extra_options):
             ):
                 raise pybamm.OptionError(
                     "If there are multiple particle phases: 'surface form' cannot be "
-                    "'false', 'particle size' must be 'false', 'particle' must be "
+                    "'false', 'particle size' must be 'single', 'particle' must be "
                     "'Fickian diffusion'. Also the following must "
                     "be 'none': 'particle mechanics', "
                     "'loss of active material', 'lithium plating'"
@@ -554,9 +555,10 @@ def __init__(self, extra_options):
                                 "interface utilisation",
                                 "loss of active material",
                                 "open circuit potential",
-                                "particle mechanics",
                                 "particle",
+                                "particle mechanics",
                                 "particle phases",
+                                "particle size",
                                 "stress-induced diffusion",
                             ]
                             and isinstance(value, tuple)

pybamm/models/submodels/active_material/base_active_material.py

@@ -83,10 +83,11 @@ def _get_standard_active_material_variables(self, eps_solid):
             # R_n, R_p. For a size distribution, calculate the area-weighted
             # mean using the distribution instead. Then the surface area is
             # calculated the same way
-            if self.options["particle size"] == "single":
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "single":
                 R = self.phase_param.R
                 R_dim = self.phase_param.R_dimensional
-            elif self.options["particle size"] == "distribution":
+            elif domain_options["particle size"] == "distribution":
                 if self.domain == "negative":
                     R_ = pybamm.standard_spatial_vars.R_n
                 elif self.domain == "positive":

pybamm/models/submodels/interface/base_interface.py

@@ -73,7 +73,8 @@ def _get_exchange_current_density(self, variables):
         if self.reaction == "lithium-ion main":
             # For "particle-size distribution" submodels, take distribution version
             # of c_s_surf that depends on particle size.
-            if self.options["particle size"] == "distribution":
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
                 c_s_surf = variables[
                     f"{Domain} {phase_name}particle surface concentration distribution"
                 ]

pybamm/models/submodels/interface/kinetics/base_kinetics.py

@@ -72,16 +72,17 @@ def get_coupled_variables(self, variables):
                 delta_phi = delta_phi.orphans[0]
         # For "particle-size distribution" models, delta_phi must then be
         # broadcast to "particle size" domain
+        domain_options = getattr(self.options, domain)
         if (
             self.reaction == "lithium-ion main"
-            and self.options["particle size"] == "distribution"
+            and domain_options["particle size"] == "distribution"
         ):
             delta_phi = pybamm.PrimaryBroadcast(delta_phi, [f"{domain} particle size"])
 
         # Get exchange-current density
         j0 = self._get_exchange_current_density(variables)
         # Get open-circuit potential variables and reaction overpotential
-        if self.options["particle size"] == "distribution":
+        if domain_options["particle size"] == "distribution":
             ocp = variables[
                 f"{Domain} electrode {reaction_name}open circuit potential distribution"
             ]

pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py

@@ -12,13 +12,26 @@ def get_coupled_variables(self, variables):
         m_lith = pybamm.sigmoid(current, 0, k)  # for lithation (current < 0)
         m_delith = 1 - m_lith  # for delithiation (current > 0)
 
-        Domain = self.domain.capitalize()
+        domain, Domain = self.domain_Domain
         phase_name = self.phase_name
 
         if self.reaction == "lithium-ion main":
             T = variables[f"{Domain} electrode temperature"]
-            # Particle size distribution is not yet implemented
-            if self.options["particle size"] != "distribution":
+            # For "particle-size distribution" models, take distribution version
+            # of c_s_surf that depends on particle size.
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
+                c_s_surf = variables[
+                    f"{Domain} {phase_name}particle surface concentration distribution"
+                ]
+                # If variable was broadcast, take only the orphan
+                if isinstance(c_s_surf, pybamm.Broadcast) and isinstance(
+                    T, pybamm.Broadcast
+                ):
+                    c_s_surf = c_s_surf.orphans[0]
+                    T = T.orphans[0]
+                T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"])
+            else:
                 c_s_surf = variables[
                     f"{Domain} {phase_name}particle surface concentration"
                 ]

pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py

@@ -11,11 +11,11 @@ def get_coupled_variables(self, variables):
         phase_name = self.phase_name
 
         if self.reaction == "lithium-ion main":
-
             T = variables[f"{Domain} electrode temperature"]
             # For "particle-size distribution" models, take distribution version
             # of c_s_surf that depends on particle size.
-            if self.options["particle size"] == "distribution":
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
                 c_s_surf = variables[
                     f"{Domain} {phase_name}particle surface concentration distribution"
                 ]

pybamm/models/submodels/particle/base_particle.py

@@ -26,7 +26,8 @@ class BaseParticle(pybamm.BaseSubModel):
     def __init__(self, param, domain, options, phase="primary"):
         super().__init__(param, domain, options=options, phase=phase)
         # Read from options to see if we have a particle size distribution
-        self.size_distribution = self.options["particle size"] == "distribution"
+        domain_options = getattr(self.options, domain)
+        self.size_distribution = domain_options["particle size"] == "distribution"
 
     def _get_effective_diffusivity(self, c, T):
         param = self.param

pybamm/parameters/size_distribution_parameters.py

@@ -17,6 +17,7 @@ def get_size_distribution_parameters(
     R_min_p=None,
     R_max_n=None,
     R_max_p=None,
+    electrode="both",
 ):
     """
     A convenience method to add standard area-weighted particle-size distribution
@@ -54,50 +55,70 @@ def get_size_distribution_parameters(
     R_max_p : float (optional)
         Maximum radius in positive electrode, scaled by the mean radius R_p_av.
         Default is 5 standard deviations above the mean.
+    electrode : str (optional)
+        Which electrode to add parameters for. If "both" (default), size distribution
+        parameters are added for both electrodes. Otherwise can be "negative" or
+        "positive" to indicate a half-cell model, in which case size distribution
+        parameters are only added for a single electrode.
     """
-
-    # Radii from given parameter set
-    R_n_typ = param["Negative particle radius [m]"]
-    R_p_typ = param["Positive particle radius [m]"]
-
-    # Set the mean particle radii for each electrode
-    R_n_av = R_n_av or R_n_typ
-    R_p_av = R_p_av or R_p_typ
-
-    # Minimum radii
-    R_min_n = R_min_n or np.max([0, 1 - sd_n * 5])
-    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)
-    R_max_p = R_max_p or (1 + sd_p * 5)
-
-    # Area-weighted particle-size distributions
-    def f_a_dist_n_dim(R):
-        return lognormal(R, R_n_av, sd_n * R_n_av)
-
-    def f_a_dist_p_dim(R):
-        return lognormal(R, R_p_av, sd_p * R_p_av)
-
-    param.update(
-        {
-            "Negative area-weighted mean particle radius [m]": R_n_av,
-            "Positive area-weighted mean particle radius [m]": R_p_av,
-            "Negative area-weighted particle-size "
-            + "standard deviation [m]": sd_n * R_n_av,
-            "Positive area-weighted particle-size "
-            + "standard deviation [m]": sd_p * R_p_av,
-            "Negative minimum particle radius [m]": R_min_n * R_n_av,
-            "Positive minimum particle radius [m]": R_min_p * R_p_av,
-            "Negative maximum particle radius [m]": R_max_n * R_n_av,
-            "Positive maximum particle radius [m]": R_max_p * R_p_av,
-            "Negative area-weighted "
-            + "particle-size distribution [m-1]": f_a_dist_n_dim,
-            "Positive area-weighted "
-            + "particle-size distribution [m-1]": f_a_dist_p_dim,
-        },
-        check_already_exists=False,
-    )
+    if electrode in ["both", "negative"]:
+        # 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_dim(R):
+            return lognormal(R, R_n_av, sd_n * R_n_av)
+
+        param.update(
+            {
+                "Negative area-weighted mean particle radius [m]": R_n_av,
+                "Negative area-weighted particle-size "
+                + "standard deviation [m]": sd_n * R_n_av,
+                "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_dim,
+            },
+            check_already_exists=False,
+        )
+    if electrode in ["both", "positive"]:
+        # 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_dim(R):
+            return lognormal(R, R_p_av, sd_p * R_p_av)
+
+        param.update(
+            {
+                "Positive area-weighted mean particle radius [m]": R_p_av,
+                "Positive area-weighted particle-size "
+                + "standard deviation [m]": sd_p * R_p_av,
+                "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_dim,
+            },
+            check_already_exists=False,
+        )
     return param
 
 

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

@@ -59,6 +59,14 @@ def test_basic_processing(self):
         )
         modeltest.test_all()
 
+    def test_basic_processing_tuple(self):
+        options = {"particle size": ("single", "distribution")}
+        model = pybamm.lithium_ion.DFN(options)
+        modeltest = tests.StandardModelTest(
+            model, parameter_values=self.params, var_pts=self.var_pts
+        )
+        modeltest.test_all()
+
     def test_uniform_profile(self):
         options = {"particle size": "distribution", "particle": "uniform profile"}
         model = pybamm.lithium_ion.DFN(options)

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

@@ -23,6 +23,10 @@ def test_well_posed_size_distribution_uniform_profile(self):
         options = {"particle size": "distribution", "particle": "uniform profile"}
         self.check_well_posedness(options)
 
+    def test_well_posed_size_distribution_tuple(self):
+        options = {"particle size": ("single", "distribution")}
+        self.check_well_posedness(options)
+
     def test_well_posed_external_circuit_explicit_power(self):
         options = {"operating mode": "explicit power"}
         self.check_well_posedness(options)

tests/unit/test_parameters/test_size_distribution_parameters.py

@@ -13,8 +13,15 @@ def test_parameter_values(self):
         values = pybamm.lithium_ion.BaseModel().default_parameter_values
         param = pybamm.LithiumIonParameters()
 
-        # add distribution parameter values
-        values = pybamm.get_size_distribution_parameters(values)
+        # add distribution parameter values for negative electrode
+        values = pybamm.get_size_distribution_parameters(values, electrode="negative")
+
+        # check positive parameters aren't there yet 
+        with self.assertRaises(KeyError):
+            values["Positive maximum particle radius [m]"]
+
+        # now add distribution parameter values for positive electrode
+        values = pybamm.get_size_distribution_parameters(values, electrode="positive")
 
         # check dimensionless parameters
 
@@ -34,6 +41,7 @@ def test_parameter_values(self):
         values.evaluate(param.p.prim.f_a_dist(R_test))
 
 
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys