Issue 3549 psd diffusivity bug (#4726)

* add tertiary integral to FV method * remove unreachable line * #3549 update flux calcs * #3549 fix averaging, add test * #3549 changelog * #3549 coverage * #3549 fix N_s_xav in x-averaged PSD case --------- Co-authored-by: Eric G. Kratz <[email protected]> Co-authored-by: Alexander Bills <[email protected]> Co-authored-by: Alec Bills <[email protected]>
pybamm-team · Jan 6, 2025 · b0a5c31 · b0a5c31
1 parent 8ec98fe
commit b0a5c31
Show file tree

Hide file tree

Showing 6 changed files with 56 additions and 29 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -22,6 +22,7 @@ package to install PyBaMM with only the required dependencies. ([conda-forge/pyb
 
 ## Bug fixes
 
+- Fixed bug when using stoichiometry-dependent diffusivity with the DFN model with a particle size distribution. ([#4726](https://github.com/pybamm-team/PyBaMM/pull/4726))
 - Remove internal use of deprecated `set_parameters` function in the `Simulation` class which caused warnings. ([#4638](https://github.com/pybamm-team/PyBaMM/pull/4638))
 - Provide default value for `Symbol.mesh` attribute to avoid errors when adding variables after discretisation. ([#4644](https://github.com/pybamm-team/PyBaMM/pull/4644))
 

diff --git a/src/pybamm/models/submodels/particle/base_particle.py b/src/pybamm/models/submodels/particle/base_particle.py
@@ -166,7 +166,7 @@ def _get_standard_flux_variables(self, N_s):
 
         variables = {f"{Domain} {phase_name}particle flux [mol.m-2.s-1]": N_s}
 
-        if isinstance(N_s, pybamm.Broadcast):
+        if isinstance(N_s, pybamm.SecondaryBroadcast):
             N_s_xav = pybamm.x_average(N_s)
             variables.update(
                 {
@@ -403,29 +403,29 @@ def _get_standard_flux_distribution_variables(self, N_s):
 
         if [f"{domain} electrode"] in N_s.domains.values():
             # N_s depends on x
-
             N_s_distribution = N_s
-            # x-av the *tertiary* domain
-            # NOTE: not yet implemented. Fill with zeros instead
-            N_s_xav_distribution = pybamm.FullBroadcast(
-                0,
-                [f"{domain} {phase_name}particle"],
-                {
-                    "secondary": f"{domain} {phase_name}particle size",
-                    "tertiary": "current collector",
-                },
-            )
+            if isinstance(N_s, pybamm.TertiaryBroadcast):
+                N_s_xav_distribution = N_s.orphans[0]
+            else:
+                # can't average variables that evaluate on edges
+                N_s_xav_distribution = None
         else:
             N_s_xav_distribution = N_s
             N_s_distribution = pybamm.TertiaryBroadcast(N_s, [f"{domain} electrode"])
 
         variables = {
-            f"X-averaged {domain} {phase_name}particle flux "
-            "distribution [mol.m-2.s-1]": N_s_xav_distribution,
             f"{Domain} {phase_name}particle flux "
             "distribution [mol.m-2.s-1]": N_s_distribution,
         }
 
+        if N_s_xav_distribution is not None:
+            variables.update(
+                {
+                    f"X-averaged {domain} {phase_name}particle flux "
+                    "distribution [mol.m-2.s-1]": N_s_xav_distribution,
+                }
+            )
+
         return variables
 
     def _get_standard_diffusivity_variables(self, D_eff):

diff --git a/src/pybamm/models/submodels/particle/fickian_diffusion.py b/src/pybamm/models/submodels/particle/fickian_diffusion.py
@@ -130,7 +130,6 @@ 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:
@@ -213,22 +212,25 @@ def get_coupled_variables(self, variables):
             }
         )
 
+        if self.x_average is True:
+            if self.size_distribution is True:
+                D_eff = pybamm.TertiaryBroadcast(D_eff, [f"{domain} electrode"])
+                N_s = pybamm.TertiaryBroadcast(N_s, [f"{domain} electrode"])
+            else:
+                D_eff = pybamm.SecondaryBroadcast(D_eff, [f"{domain} electrode"])
+                N_s = pybamm.SecondaryBroadcast(N_s, [f"{domain} electrode"])
+
         if self.size_distribution is True:
             # Size-dependent flux variables
             variables.update(
                 self._get_standard_diffusivity_distribution_variables(D_eff)
             )
             variables.update(self._get_standard_flux_distribution_variables(N_s))
             # Size-averaged flux variables
-            R = variables[f"{Phase_prefactor}{Domain} {phase_name}particle sizes [m]"]
+            D_eff = pybamm.size_average(D_eff)
+            R = pybamm.SpatialVariable("R", domains=D_eff.domains)
             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)
-
-        if self.x_average is True:
-            D_eff = pybamm.SecondaryBroadcast(D_eff, [f"{domain} electrode"])
-            N_s = pybamm.SecondaryBroadcast(N_s, [f"{domain} electrode"])
-
         variables.update(self._get_standard_diffusivity_variables(D_eff))
         variables.update(self._get_standard_flux_variables(N_s))
 

diff --git a/src/pybamm/models/submodels/particle/msmr_diffusion.py b/src/pybamm/models/submodels/particle/msmr_diffusion.py
@@ -240,22 +240,25 @@ def get_coupled_variables(self, variables):
             }
         )
 
+        if self.x_average is True:
+            if self.size_distribution is True:
+                D_eff = pybamm.TertiaryBroadcast(D_eff, [f"{domain} electrode"])
+                N_s = pybamm.TertiaryBroadcast(N_s, [f"{domain} electrode"])
+            else:
+                D_eff = pybamm.SecondaryBroadcast(D_eff, [f"{domain} electrode"])
+                N_s = pybamm.SecondaryBroadcast(N_s, [f"{domain} electrode"])
+
         if self.size_distribution is True:
             # Size-dependent flux variables
             variables.update(
                 self._get_standard_diffusivity_distribution_variables(D_eff)
             )
             variables.update(self._get_standard_flux_distribution_variables(N_s))
             # Size-averaged flux variables
-            R = variables[f"{Domain} {phase_name}particle sizes [m]"]
+            D_eff = pybamm.size_average(D_eff)
+            R = pybamm.SpatialVariable("R", domains=D_eff.domains)
             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)
-
-        if self.x_average is True:
-            D_eff = pybamm.SecondaryBroadcast(D_eff, [f"{domain} electrode"])
-            N_s = pybamm.SecondaryBroadcast(N_s, [f"{domain} electrode"])
-
         variables.update(self._get_standard_diffusivity_variables(D_eff))
         variables.update(self._get_standard_flux_variables(N_s))
 

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
@@ -151,3 +151,16 @@ def test_conservation_each_electrode(self):
         # compare
         np.testing.assert_array_almost_equal(neg_Li[0], neg_Li[1], decimal=12)
         np.testing.assert_array_almost_equal(pos_Li[0], pos_Li[1], decimal=12)
+
+    def test_basic_processing_nonlinear_diffusion(self):
+        options = {
+            "particle size": "distribution",
+        }
+        model = pybamm.lithium_ion.DFN(options)
+        # Ecker2015 has a nonlinear diffusion coefficient
+        parameter_values = pybamm.ParameterValues("Ecker2015")
+        parameter_values = pybamm.get_size_distribution_parameters(parameter_values)
+        modeltest = tests.StandardModelTest(
+            model, parameter_values=parameter_values, var_pts=self.var_pts
+        )
+        modeltest.test_all(skip_output_tests=True)
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
@@ -122,3 +122,11 @@ def test_conservation_each_electrode(self):
         # compare
         np.testing.assert_array_almost_equal(neg_Li[0], neg_Li[1], decimal=13)
         np.testing.assert_array_almost_equal(pos_Li[0], pos_Li[1], decimal=13)
+
+    def test_basic_processing_nonlinear_diffusion(self):
+        model = pybamm.lithium_ion.MPM()
+        # Ecker2015 has a nonlinear diffusion coefficient
+        parameter_values = pybamm.ParameterValues("Ecker2015")
+        parameter_values = pybamm.get_size_distribution_parameters(parameter_values)
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all(skip_output_tests=True)