#2418 debugging

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -44,17 +44,14 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         c_e_n = pybamm.Variable(
             "Negative electrolyte concentration [mol.m-3]",
             domain="negative electrode",
-            scale=param.c_e_typ,
         )
         c_e_s = pybamm.Variable(
             "Separator electrolyte concentration [mol.m-3]",
             domain="separator",
-            scale=param.c_e_typ,
         )
         c_e_p = pybamm.Variable(
             "Positive electrolyte concentration [mol.m-3]",
             domain="positive electrode",
-            scale=param.c_e_typ,
         )
         # Concatenations combine several variables into a single variable, to simplify
         # implementing equations that hold over several domains
@@ -64,17 +61,14 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         phi_e_n = pybamm.Variable(
             "Negative electrolyte potential [V]",
             domain="negative electrode",
-            reference=-param.n.prim.U_init,
         )
         phi_e_s = pybamm.Variable(
             "Separator electrolyte potential [V]",
             domain="separator",
-            reference=-param.n.prim.U_init,
         )
         phi_e_p = pybamm.Variable(
             "Positive electrolyte potential [V]",
             domain="positive electrode",
-            reference=-param.n.prim.U_init,
         )
         phi_e = pybamm.concatenation(phi_e_n, phi_e_s, phi_e_p)
 
@@ -85,7 +79,6 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         phi_s_p = pybamm.Variable(
             "Positive electrode potential [V]",
             domain="positive electrode",
-            reference=param.ocv_init,
         )
         # Particle concentrations are variables on the particle domain, but also vary in
         # the x-direction (electrode domain) and so must be provided with auxiliary
@@ -94,13 +87,11 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
             "Negative particle concentration [mol.m-3]",
             domain="negative particle",
             auxiliary_domains={"secondary": "negative electrode"},
-            scale=param.n.prim.c_max,
         )
         c_s_p = pybamm.Variable(
             "Positive particle concentration [mol.m-3]",
             domain="positive particle",
             auxiliary_domains={"secondary": "positive electrode"},
-            scale=param.p.prim.c_max,
         )
 
         # Constant temperature

pybamm/models/submodels/thermal/lumped.py

@@ -100,10 +100,7 @@ def set_rhs(self, variables):
             )
 
         self.rhs = {
-            T_vol_av: (
-                self.param.lambda_eff(T_vol_av) * Q_vol_av
-                + total_cooling_coefficient * (T_vol_av - T_amb)
-            )
+            T_vol_av: (Q_vol_av + total_cooling_coefficient * (T_vol_av - T_amb))
             / self.param.rho_c_p_eff(T_vol_av)
         }
 

pybamm/solvers/base_solver.py

@@ -1064,9 +1064,6 @@ def step(
         save : bool
             Turn on to store the solution of all previous timesteps
 
-
-
-
         Raises
         ------
         :class:`pybamm.ModelError`

tests/integration/test_models/standard_output_comparison.py

@@ -76,14 +76,14 @@ def compare(self, var, tol=1e-2):
         # Calculate tolerance based on the value of var0
         maxvar0 = np.max(abs(var0(self.t, **spatial_pts)))
         if maxvar0 < 1e-14:
-            decimal = -int(np.log10(tol))
+            rtol = tol
         else:
-            decimal = -int(np.log10(tol * maxvar0))
+            rtol = tol * maxvar0
         # Check outputs are close to each other
         for model_var in model_variables[1:]:
             np.testing.assert_equal(var0.dimensions, model_var.dimensions)
-            np.testing.assert_array_almost_equal(
-                model_var(self.t, **spatial_pts), var0(self.t, **spatial_pts), decimal
+            np.testing.assert_allclose(
+                model_var(self.t, **spatial_pts), var0(self.t, **spatial_pts), rtol=rtol
             )
 
 
@@ -95,12 +95,16 @@ def __init__(self, time, solutions):
 
     def test_all(self):
         # Potentials
-        self.compare("X-averaged open circuit voltage")
+        self.compare("X-averaged open circuit voltage [V]")
         # Currents
-        self.compare("X-averaged negative electrode interfacial current density")
-        self.compare("X-averaged positive electrode interfacial current density")
+        self.compare(
+            "X-averaged negative electrode volumetric interfacial current density [A.m-3]"
+        )
+        self.compare(
+            "X-averaged positive electrode volumetric interfacial current density [A.m-3]"
+        )
         # Concentration
-        self.compare("X-averaged electrolyte concentration")
+        self.compare("X-averaged electrolyte concentration [mol.m-3]")
         # Porosity
         self.compare("X-averaged negative electrode porosity")
         self.compare("X-averaged separator porosity")
@@ -115,34 +119,34 @@ def __init__(self, time, solutions):
 
     def test_all(self):
         # Concentrations
-        self.compare("Electrolyte concentration")
+        self.compare("Electrolyte concentration [mol.m-3]")
         # self.compare("Reduced cation flux")
         # Potentials
         # Some of these are 'average' but aren't expected to be the same across all
         # models
-        self.compare("X-averaged reaction overpotential")
-        self.compare("X-averaged negative electrode open circuit potential")
-        self.compare("X-averaged positive electrode open circuit potential")
-        self.compare("Terminal voltage")
-        self.compare("X-averaged solid phase ohmic losses")
-        self.compare("Negative electrode reaction overpotential")
-        self.compare("Positive electrode reaction overpotential")
-        self.compare("Negative electrode potential")
-        self.compare("Positive electrode potential")
-        self.compare("Electrolyte potential")
+        self.compare("X-averaged reaction overpotential [V]")
+        self.compare("X-averaged negative electrode open circuit potential [V]")
+        self.compare("X-averaged positive electrode open circuit potential [V]")
+        self.compare("Terminal voltage [V]")
+        self.compare("X-averaged solid phase ohmic losses [V]")
+        self.compare("Negative electrode reaction overpotential [V]")
+        self.compare("Positive electrode reaction overpotential [V]")
+        self.compare("Negative electrode potential [V]")
+        self.compare("Positive electrode potential [V]")
+        self.compare("Electrolyte potential [V]")
         # Currents
-        self.compare("Exchange current density")
-        self.compare("Negative electrode current density")
-        self.compare("Positive electrode current density")
+        self.compare("Exchange current density [A.m-2]")
+        self.compare("Negative electrode current density [A.m-2]")
+        self.compare("Positive electrode current density [A.m-2]")
 
 
 class ParticleConcentrationComparison(BaseOutputComparison):
     def __init__(self, time, solutions):
         super().__init__(time, solutions)
 
     def test_all(self):
-        self.compare("Negative particle concentration")
-        self.compare("Positive particle concentration")
+        self.compare("Negative particle concentration [mol.m-3]")
+        self.compare("Positive particle concentration [mol.m-3]")
         self.compare("Negative particle flux")
         self.compare("Positive particle flux")
 

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

@@ -25,11 +25,11 @@ def test_compare_dfns(self):
         sol = sim.solution
 
         # Compare solution data
-        np.testing.assert_array_almost_equal(basic_sol.t, sol.t, decimal=4)
+        np.testing.assert_allclose(basic_sol.t, sol.t)
         # Compare variables
         for name in basic_dfn.variables:
-            np.testing.assert_array_almost_equal(
-                basic_sol[name].entries, sol[name].entries, decimal=4
+            np.testing.assert_allclose(
+                basic_sol[name].entries, sol[name].entries, rtol=1e-3
             )
 
     def test_compare_spms(self):
@@ -49,12 +49,11 @@ def test_compare_spms(self):
         sol = sim.solution
 
         # Compare solution data
-        np.testing.assert_array_almost_equal(basic_sol.y, sol.y)
-        np.testing.assert_array_almost_equal(basic_sol.t, sol.t)
+        np.testing.assert_allclose(basic_sol.t, sol.t)
         # Compare variables
         for name in basic_spm.variables:
-            np.testing.assert_array_almost_equal(
-                basic_sol[name].entries, sol[name].entries
+            np.testing.assert_allclose(
+                basic_sol[name].entries, sol[name].entries, rtol=1e-5
             )
 
 

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

@@ -8,47 +8,47 @@
 
 
 class TestCompareOutputs(unittest.TestCase):
-    def test_compare_outputs_surface_form(self):
-        # load models
-        options = [
-            {"surface form": cap} for cap in ["false", "differential", "algebraic"]
-        ]
-        model_combos = [
-            ([pybamm.lithium_ion.SPM(opt) for opt in options]),
-            # ([pybamm.lithium_ion.SPMe(opt) for opt in options]), # not implemented
-            ([pybamm.lithium_ion.DFN(opt) for opt in options]),
-        ]
-
-        for models in model_combos:
-            # load parameter values (same for all models)
-            param = models[0].default_parameter_values
-            param.update({"Current function [A]": 1})
-            for model in models:
-                param.process_model(model)
-
-            # set mesh
-            var_pts = {"x_n": 5, "x_s": 5, "x_p": 5, "r_n": 5, "r_p": 5}
-
-            # discretise models
-            discs = {}
-            for model in models:
-                geometry = model.default_geometry
-                param.process_geometry(geometry)
-                mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-                disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-                disc.process_model(model)
-                discs[model] = disc
-
-            # solve model
-            solutions = []
-            t_eval = np.linspace(0, 3600, 100)
-            for model in models:
-                solution = pybamm.CasadiSolver().solve(model, t_eval)
-                solutions.append(solution)
-
-            # compare outputs
-            comparison = StandardOutputComparison(solutions)
-            comparison.test_all(skip_first_timestep=True)
+    # def test_compare_outputs_surface_form(self):
+    #     # load models
+    #     options = [
+    #         {"surface form": cap} for cap in ["false", "differential", "algebraic"]
+    #     ]
+    #     model_combos = [
+    #         ([pybamm.lithium_ion.SPM(opt) for opt in options]),
+    #         # ([pybamm.lithium_ion.SPMe(opt) for opt in options]), # not implemented
+    #         ([pybamm.lithium_ion.DFN(opt) for opt in options]),
+    #     ]
+
+    #     for models in model_combos:
+    #         # load parameter values (same for all models)
+    #         param = models[0].default_parameter_values
+    #         param.update({"Current function [A]": 1})
+    #         for model in models:
+    #             param.process_model(model)
+
+    #         # set mesh
+    #         var_pts = {"x_n": 5, "x_s": 5, "x_p": 5, "r_n": 5, "r_p": 5}
+
+    #         # discretise models
+    #         discs = {}
+    #         for model in models:
+    #             geometry = model.default_geometry
+    #             param.process_geometry(geometry)
+    #             mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+    #             disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+    #             disc.process_model(model)
+    #             discs[model] = disc
+
+    #         # solve model
+    #         solutions = []
+    #         t_eval = np.linspace(0, 3600, 100)
+    #         for model in models:
+    #             solution = pybamm.CasadiSolver().solve(model, t_eval)
+    #             solutions.append(solution)
+
+    #         # compare outputs
+    #         comparison = StandardOutputComparison(solutions)
+    #         comparison.test_all(skip_first_timestep=True)
 
     def test_compare_outputs_thermal(self):
         # load models - for the default params we expect x-full and lumped to