Remove harmonic mean regularization (#4977)

* smaller magic number in `FunctionParameter`

* remove regularization in harmonic mean

change the functional form of the harmonic mean to reduce number of diffusion variables in codegen

* fix tests

* style

* Update CHANGELOG.md

Author: MarcBerliner

CHANGELOG.md

@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Bug fixes
+
+- Remove a regularization term in the harmonic mean.  ([#4977](https://github.com/pybamm-team/PyBaMM/pull/4977))
+
 # [v25.4.0](https://github.com/pybamm-team/PyBaMM/tree/v25.4.0) - 2025-04-02
 
 ## Features

src/pybamm/expression_tree/parameter.py

@@ -217,8 +217,8 @@ def _evaluate_for_shape(self):
         Returns the sum of the evaluated children
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
         """
-        # add 1e-16 to avoid division by zero
-        return sum(child.evaluate_for_shape() for child in self.children) + 1e-16
+        # add 1e-300 to avoid division by zero
+        return sum(child.evaluate_for_shape() for child in self.children) + 1e-300
 
     def to_equation(self) -> sympy.Symbol:
         """Convert the node and its subtree into a SymPy equation."""

src/pybamm/spatial_methods/finite_volume.py

@@ -1334,7 +1334,7 @@ def harmonic_mean(array):
             The harmonic mean is computed as
 
             .. math::
-                D_{eff} = \\frac{D_1  D_2}{\\beta D_2 + (1 - \\beta) D_1},
+                D_{eff} = \\frac{1}{\\frac{\\beta}{D_1} + \\frac{1 - \\beta}{D_2}},
 
             where
 
@@ -1409,8 +1409,7 @@ def harmonic_mean(array):
 
                 # dx_real = dx * length, therefore, beta is unchanged
                 # Compute harmonic mean on internal edges
-                # Note: add small number to denominator to regularise D_eff
-                D_eff = D1 * D2 / (D2 * beta + D1 * (1 - beta) + 1e-16)
+                D_eff = 1 / (beta / D1 + (1 - beta) / D2)
 
                 # Matrix to pad zeros at the beginning and end of the array where
                 # the exterior edge values will be added
@@ -1452,8 +1451,7 @@ def harmonic_mean(array):
                 beta = pybamm.Array(np.kron(np.ones((second_dim_repeats, 1)), sub_beta))
 
                 # Compute harmonic mean on nodes
-                # Note: add small number to denominator to regularise D_eff
-                D_eff = D1 * D2 / (D2 * beta + D1 * (1 - beta) + 1e-16)
+                D_eff = 1 / (beta / D1 + (1 - beta) / D2)
 
                 return D_eff
 

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

@@ -13,14 +13,17 @@ def test_compare_dfns(self):
         dfn = pybamm.lithium_ion.DFN()
 
         # Solve basic DFN
-        basic_sim = pybamm.Simulation(basic_dfn, parameter_values=parameter_values)
+        solver = pybamm.IDAKLUSolver(atol=1e-8, rtol=1e-8)
+        basic_sim = pybamm.Simulation(
+            basic_dfn, parameter_values=parameter_values, solver=solver
+        )
         t_eval = [0, 3600]
         t_interp = np.linspace(0, 3600)
         basic_sim.solve(t_eval=t_eval, t_interp=t_interp)
         basic_sol = basic_sim.solution
 
         # Solve main DFN
-        sim = pybamm.Simulation(dfn, parameter_values=parameter_values)
+        sim = pybamm.Simulation(dfn, parameter_values=parameter_values, solver=solver)
         sim.solve(t_eval=t_eval, t_interp=t_interp)
         sol = sim.solution
 
@@ -43,7 +46,10 @@ def test_compare_dfns_composite(self):
         parameter_values = pybamm.ParameterValues("Chen2020_composite")
 
         # Solve basic DFN
-        basic_sim = pybamm.Simulation(basic_dfn, parameter_values=parameter_values)
+        solver = pybamm.IDAKLUSolver(atol=1e-8, rtol=1e-8)
+        basic_sim = pybamm.Simulation(
+            basic_dfn, parameter_values=parameter_values, solver=solver
+        )
         t_eval = [0, 3600]
         t_interp = np.linspace(0, 3600)
         basic_sim.solve(t_eval=t_eval, t_interp=t_interp)
@@ -69,14 +75,17 @@ def test_compare_spms(self):
         spm = pybamm.lithium_ion.SPM()
 
         # Solve basic SPM
-        basic_sim = pybamm.Simulation(basic_spm, parameter_values=parameter_values)
+        solver = pybamm.IDAKLUSolver(atol=1e-8, rtol=1e-8)
+        basic_sim = pybamm.Simulation(
+            basic_spm, parameter_values=parameter_values, solver=solver
+        )
         t_eval = [0, 3600]
         t_interp = np.linspace(0, 3600)
         basic_sim.solve(t_eval=t_eval, t_interp=t_interp)
         basic_sol = basic_sim.solution
 
         # Solve main SPM
-        sim = pybamm.Simulation(spm, parameter_values=parameter_values)
+        sim = pybamm.Simulation(spm, parameter_values=parameter_values, solver=solver)
         sim.solve(t_eval=t_eval, t_interp=t_interp)
         sol = sim.solution
 

tests/unit/test_experiments/test_experiment.py

@@ -276,7 +276,8 @@ def test_voltage_without_directions(self):
         step = pybamm.step.voltage(2.5, termination="2.5 V")
         experiment = pybamm.Experiment([step])
 
-        sim = pybamm.Simulation(model, experiment=experiment)
+        solver = pybamm.IDAKLUSolver(atol=1e-8, rtol=1e-8)
+        sim = pybamm.Simulation(model, experiment=experiment, solver=solver)
 
         sim.solve()
         solution = sim.solution

tests/unit/test_experiments/test_simulation_with_experiment.py

@@ -73,7 +73,8 @@ def test_run_experiment(self):
             temperature="-14oC",
         )
         model = pybamm.lithium_ion.SPM()
-        sim = pybamm.Simulation(model, experiment=experiment)
+        solver = pybamm.IDAKLUSolver(atol=1e-8, rtol=1e-8)
+        sim = pybamm.Simulation(model, experiment=experiment, solver=solver)
         # test the callback here
         sol = sim.solve(callbacks=pybamm.callbacks.Callback())
         assert sol.termination == "final time"