#2418 fix unit tests, add scaling to some algebraic equations

examples/scripts/compare_lithium_ion.py

@@ -3,12 +3,12 @@
 #
 import pybamm
 
-pybamm.set_logging_level("INFO")
+pybamm.set_logging_level("DEBUG")
 # load models
 models = [
     # pybamm.lithium_ion.SPM(),
     # pybamm.lithium_ion.SPMe(),
-    pybamm.lithium_ion.DFN({"dimensionality": 2}),
+    pybamm.lithium_ion.BasicDFN(),
     # pybamm.lithium_ion.DFN(
     #     {"particle": "uniform profile"}
     # ),
@@ -20,9 +20,9 @@
 # create and run simulations
 sims = []
 for model in models:
-    sim = pybamm.Simulation(model)
-    sim.solve([0, 3600])
+    sim = pybamm.Simulation(model, solver=pybamm.CasadiSolver("fast"))
+    sim.solve([0, 4000])
     sims.append(sim)
 
 # plot
-pybamm.dynamic_plot(sims)
+# pybamm.dynamic_plot(sims)

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -198,8 +198,8 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         i_s_p = -sigma_eff_p * pybamm.grad(phi_s_p)
         # The `algebraic` dictionary contains differential equations, with the key being
         # the main scalar variable of interest in the equation
-        self.algebraic[phi_s_n] = pybamm.div(i_s_n) + a_j_n
-        self.algebraic[phi_s_p] = pybamm.div(i_s_p) + a_j_p
+        self.algebraic[phi_s_n] = param.L_x**2 * (pybamm.div(i_s_n) + a_j_n)
+        self.algebraic[phi_s_p] = param.L_x**2 * (pybamm.div(i_s_p) + a_j_p)
         self.boundary_conditions[phi_s_n] = {
             "left": (pybamm.Scalar(0), "Dirichlet"),
             "right": (pybamm.Scalar(0), "Neumann"),
@@ -220,7 +220,7 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         i_e = (param.kappa_e(c_e, T) * tor) * (
             param.chiRT_over_Fc(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
         )
-        self.algebraic[phi_e] = pybamm.div(i_e) - a_j
+        self.algebraic[phi_e] = param.L_x**2 * (pybamm.div(i_e) - a_j)
         self.boundary_conditions[phi_e] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (pybamm.Scalar(0), "Neumann"),

pybamm/models/submodels/current_collector/quite_conductive_potential_pair.py

@@ -39,7 +39,6 @@ def get_fundamental_variables(self):
         i_boundary_cc = pybamm.Variable(
             "Current collector current density [A.m-2]",
             domain="current collector",
-            scale=param.Q / (param.A_cc * param.n_electrodes_parallel),
         )
 
         variables.update(self._get_standard_current_variables(i_cc, i_boundary_cc))

pybamm/models/submodels/electrode/ohm/full_ohm.py

@@ -75,6 +75,7 @@ def set_algebraic(self, variables):
             "interfacial current densities [A.m-3]"
         ]
 
+        # multiply by Lx**2 to improve conditioning
         self.algebraic[phi_s] = self.param.L_x**2 * (pybamm.div(i_s) + sum_a_j)
 
     def set_boundary_conditions(self, variables):

pybamm/models/submodels/electrolyte_conductivity/full_conductivity.py

@@ -76,6 +76,7 @@ def set_algebraic(self, variables):
         # Override print_name
         sum_a_j.print_name = "aj"
 
+        # multiply by Lx**2 to improve conditioning
         self.algebraic = {phi_e: self.param.L_x**2 * (pybamm.div(i_e) - sum_a_j)}
 
     def set_initial_conditions(self, variables):

pybamm/solvers/base_solver.py

@@ -988,7 +988,7 @@ def _get_discontinuity_start_end_indices(self, model, inputs, t_eval):
             end_indices.append(dindex + 1)
             start_indices.append(dindex + 1)
             if dtime * (1 - eps) < t_eval[dindex] < dtime * (1 + eps):
-                t_eval[dindex] += eps
+                t_eval[dindex] *= 1 + eps
                 t_eval = np.insert(t_eval, dindex, dtime * (1 - eps))
             else:
                 t_eval = np.insert(

tests/unit/test_solvers/test_algebraic_solver.py

@@ -86,6 +86,12 @@ def algebraic_eval(self, t, y, inputs):
         ):
             solver._integrate(model, np.array([0]))
 
+        solver = pybamm.AlgebraicSolver()
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Could not find acceptable solution: solver terminated"
+        ):
+            solver._integrate(model, np.array([0]))
+
     def test_with_jacobian(self):
         A = np.array([[4, 3], [1, -1]])
         b = np.array([0, 7])

tests/unit/test_solvers/test_casadi_algebraic_solver.py

@@ -72,6 +72,11 @@ def algebraic_eval(self, t, y, inputs):
             pybamm.SolverError, "Could not find acceptable solution: .../casadi"
         ):
             solver._integrate(model, np.array([0]), {})
+        solver = pybamm.CasadiAlgebraicSolver(extra_options={"error_on_fail": False})
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Could not find acceptable solution: solver terminated"
+        ):
+            solver._integrate(model, np.array([0]), {})
 
         # Model returns Nan
         class NaNModel:
@@ -198,6 +203,7 @@ def test_least_squares_fit(self):
         # Simple system: a single algebraic equation
         var = pybamm.Variable("var", domain="negative electrode")
         model = pybamm.BaseModel()
+        model._length_scales = {"negative electrode": pybamm.Scalar(1)}
         p = pybamm.InputParameter("p")
         q = pybamm.InputParameter("q")
         model.algebraic = {var: (var - p)}
@@ -237,6 +243,7 @@ def jac(x):
     def test_solve_with_symbolic_input_1D_scalar_input(self):
         var = pybamm.Variable("var", "negative electrode")
         model = pybamm.BaseModel()
+        model._length_scales = {"negative electrode": pybamm.Scalar(1)}
         param = pybamm.InputParameter("param")
         model.algebraic = {var: var + param}
         model.initial_conditions = {var: 2}
@@ -262,6 +269,7 @@ def test_solve_with_symbolic_input_1D_scalar_input(self):
     def test_solve_with_symbolic_input_1D_vector_input(self):
         var = pybamm.Variable("var", "negative electrode")
         model = pybamm.BaseModel()
+        model._length_scales = {"negative electrode": pybamm.Scalar(1)}
         param = pybamm.InputParameter("param", "negative electrode")
         model.algebraic = {var: var + param}
         model.initial_conditions = {var: 2}
@@ -320,6 +328,7 @@ def test_least_squares_fit_input_in_initial_conditions(self):
         # Simple system: a single algebraic equation
         var = pybamm.Variable("var", domain="negative electrode")
         model = pybamm.BaseModel()
+        model._length_scales = {"negative electrode": pybamm.Scalar(1)}
         p = pybamm.InputParameter("p")
         q = pybamm.InputParameter("q")
         model.algebraic = {var: (var - p)}

tests/unit/test_solvers/test_casadi_solver.py

@@ -88,10 +88,10 @@ def test_without_grid(self):
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
         np.testing.assert_array_almost_equal(
-            solution.y[0], np.exp(0.1 * t_eval), decimal=5
+            solution.y.full()[0], np.exp(0.1 * t_eval), decimal=5
         )
         np.testing.assert_array_almost_equal(
-            solution.y[1], np.ones_like(t_eval), decimal=5
+            solution.y.full()[1], np.ones_like(t_eval), decimal=5
         )
 
         # DAE model, errors
@@ -440,7 +440,6 @@ def test_model_solver_dae_inputs_in_initial_conditions(self):
             solution.y.full()[-1], 1 * np.exp(-0.1 * solution.t), decimal=5
         )
 
-    @unittest.skip("External variables will be removed")
     def test_model_solver_with_external(self):
         # Create model
         model = pybamm.BaseModel()

tests/unit/test_solvers/test_scikits_solvers.py

@@ -339,8 +339,8 @@ def nonsmooth_mult(t):
             dindex = np.searchsorted(solution.t, discontinuity)
             value_before = solution.t[dindex - 1]
             value_after = solution.t[dindex]
-            self.assertEqual(value_before + sys.float_info.epsilon, discontinuity)
-            self.assertEqual(value_after - sys.float_info.epsilon, discontinuity)
+            self.assertEqual(value_before / (1 - sys.float_info.epsilon), discontinuity)
+            self.assertEqual(value_after / (1 + sys.float_info.epsilon), discontinuity)
 
         # both solution time vectors should have same number of points
         self.assertEqual(len(solution1.t), len(solution2.t))
@@ -406,8 +406,12 @@ def test_model_solver_dae_multiple_nonsmooth_python(self):
                 dindex = np.searchsorted(solution.t, discontinuity)
                 value_before = solution.t[dindex - 1]
                 value_after = solution.t[dindex]
-                self.assertEqual(value_before + sys.float_info.epsilon, discontinuity)
-                self.assertEqual(value_after - sys.float_info.epsilon, discontinuity)
+                self.assertEqual(
+                    value_before / (1 - sys.float_info.epsilon), discontinuity
+                )
+                self.assertEqual(
+                    value_after / (1 + sys.float_info.epsilon), discontinuity
+                )
 
         # both solution time vectors should have same number of points
         self.assertEqual(len(solution1.t), len(solution2.t))