Merge branch 'master' into issue-646-install

examples/scripts/compare_comsol/compare_comsol_DFN.py

@@ -51,7 +51,7 @@
 
 # solve model at comsol times
 time = comsol_variables["time"] / tau.evaluate(0)
-solution = pybamm_model.default_solver.solve(pybamm_model, time)
+solution = pybamm.CasadiSolver(mode="fast").solve(pybamm_model, time)
 
 
 # Make Comsol 'model' for comparison
@@ -60,12 +60,13 @@
 L_x = param.evaluate(pybamm.standard_parameters_lithium_ion.L_x)
 
 
-def get_interp_fun(variable, domain):
+def get_interp_fun(variable_name, domain):
     """
     Create a :class:`pybamm.Function` object using the variable, to allow plotting with
     :class:`'pybamm.QuickPlot'` (interpolate in space to match edges, and then create
     function to interpolate in time)
     """
+    variable = comsol_variables[variable_name]
     if domain == ["negative electrode"]:
         comsol_x = comsol_variables["x_n"]
     elif domain == ["positive electrode"]:
@@ -80,17 +81,17 @@ def myinterp(t):
         return interp.interp1d(comsol_t, variable)(t)[:, np.newaxis]
 
     # Make sure to use dimensional time
-    fun = pybamm.Function(myinterp, pybamm.t * tau)
+    fun = pybamm.Function(myinterp, pybamm.t * tau, name=variable_name + "_comsol")
     fun.domain = domain
     return fun
 
 
-comsol_c_n_surf = get_interp_fun(comsol_variables["c_n_surf"], ["negative electrode"])
-comsol_c_e = get_interp_fun(comsol_variables["c_e"], whole_cell)
-comsol_c_p_surf = get_interp_fun(comsol_variables["c_p_surf"], ["positive electrode"])
-comsol_phi_n = get_interp_fun(comsol_variables["phi_n"], ["negative electrode"])
-comsol_phi_e = get_interp_fun(comsol_variables["phi_e"], whole_cell)
-comsol_phi_p = get_interp_fun(comsol_variables["phi_p"], ["positive electrode"])
+comsol_c_n_surf = get_interp_fun("c_n_surf", ["negative electrode"])
+comsol_c_e = get_interp_fun("c_e", whole_cell)
+comsol_c_p_surf = get_interp_fun("c_p_surf", ["positive electrode"])
+comsol_phi_n = get_interp_fun("phi_n", ["negative electrode"])
+comsol_phi_e = get_interp_fun("phi_e", whole_cell)
+comsol_phi_p = get_interp_fun("phi_p", ["positive electrode"])
 comsol_voltage = interp.interp1d(comsol_t, comsol_variables["voltage"])
 
 # Create comsol model with dictionary of Matrix variables

pybamm/expression_tree/interpolant.py

@@ -52,7 +52,7 @@ def __init__(
         else:
             raise ValueError("interpolator '{}' not recognised".format(interpolator))
         # Set name
-        if name is not None:
+        if name is not None and not name.startswith("interpolating function"):
             name = "interpolating function ({})".format(name)
         else:
             name = "interpolating function"

pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py

@@ -48,7 +48,7 @@ def default_solver(self):
             return pybamm.ScipySolver()
         elif pybamm.have_scikits_odes():
             return pybamm.ScikitsDaeSolver()
-        else:
+        else:  # pragma: no cover
             return pybamm.CasadiSolver(mode="safe")
 
     def set_standard_output_variables(self):

pybamm/parameters/standard_current_functions/user_current.py

@@ -25,7 +25,7 @@ def __init__(self, function, **kwargs):
         self.function = function
 
     def __str__(self):
-        return "User defined current"
+        return "User defined current ({})".format(self.function.__name__)
 
     def __call__(self, t):
         return self.function(t, **self.parameters_eval)

pybamm/solvers/dae_solver.py

@@ -282,10 +282,10 @@ def jacobian_alg(t, y):
         def rhs(t, y):
             return concatenated_rhs_fn(t, y).full()[:, 0]
 
-        if len(model.algebraic) > 0:
+        def algebraic(t, y):
+            return concatenated_algebraic_fn(t, y).full()[:, 0]
 
-            def algebraic(t, y):
-                return concatenated_algebraic_fn(t, y).full()[:, 0]
+        if len(model.algebraic) > 0:
 
             y0 = self.calculate_consistent_initial_conditions(
                 rhs,
@@ -294,10 +294,7 @@ def algebraic(t, y):
                 jacobian_alg,
             )
         else:
-            # can use DAE solver to solve ODE model (just return empty algebraic)
-            def algebraic(t, y):
-                return np.empty(0)
-
+            # can use DAE solver to solve ODE model
             y0 = model.concatenated_initial_conditions[:, 0]
 
         # Create functions to evaluate residuals

pybamm/solvers/idaklu_solver.py

@@ -211,7 +211,7 @@ def rootfn(t, y):
 
         # get ids of rhs and algebraic variables
         rhs_ids = np.ones(self.rhs(0, y0).shape)
-        alg_ids = np.zeros(self.algebraic(0, y0).shape)
+        alg_ids = np.zeros(len(y0) - len(rhs_ids))
         ids = np.concatenate((rhs_ids, alg_ids))
 
         # solve

tests/unit/test_expression_tree/test_interpolant.py

@@ -72,6 +72,15 @@ def test_diff(self):
                 decimal=3,
             )
 
+    def test_processing(self):
+        x = np.linspace(0, 1)[:, np.newaxis]
+        y = pybamm.StateVector(slice(0, 2))
+        linear = np.hstack([x, 2 * x])
+        interp = pybamm.Interpolant(linear, y)
+
+        self.assertEqual(interp.id, interp.new_copy().id)
+        self.assertEqual(interp.id, interp.simplify().id)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

tests/unit/test_parameters/test_current_functions.py

@@ -60,6 +60,7 @@ def my_fun(t, A, omega):
         # pass my_fun to UserCurrent class, giving the additonal parameters as
         # keyword arguments
         current = pybamm.UserCurrent(my_fun, A=A, omega=omega)
+        self.assertEqual(str(current), "User defined current (my_fun)")
 
         # set and process parameters
         parameter_values = pybamm.ParameterValues(

tests/unit/test_solvers/test_scikits_solvers.py

@@ -394,9 +394,9 @@ def jacobian(t, y):
         np.testing.assert_allclose(0.125 * solution.t, solution.y[0])
         np.testing.assert_allclose(0.25 * solution.t, solution.y[1])
 
-    def test_model_solver_ode(self):
-        # Create model
+    def test_model_solver_ode_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=whole_cell)
         model.rhs = {var: 0.1 * var}
@@ -411,9 +411,9 @@ def test_model_solver_ode(self):
         np.testing.assert_array_equal(solution.t, t_eval)
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
 
-    def test_model_solver_ode_events(self):
-        # Create model
+    def test_model_solver_ode_events_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=whole_cell)
         model.rhs = {var: 0.1 * var}
@@ -433,9 +433,9 @@ def test_model_solver_ode_events(self):
         np.testing.assert_array_less(solution.y[0], 1.5)
         np.testing.assert_array_less(solution.y[0], 1.25)
 
-    def test_model_solver_ode_jacobian(self):
-        # Create model
+    def test_model_solver_ode_jacobian_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var1 = pybamm.Variable("var1", domain=whole_cell)
         var2 = pybamm.Variable("var2", domain=whole_cell)
@@ -482,9 +482,9 @@ def jacobian(t, y):
             np.ones((N, T.size)) * (T[np.newaxis, :] - np.exp(T[np.newaxis, :])),
         )
 
-    def test_model_solver_dae(self):
-        # Create model
+    def test_model_solver_dae_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var1 = pybamm.Variable("var1", domain=whole_cell)
         var2 = pybamm.Variable("var2", domain=whole_cell)
@@ -503,9 +503,9 @@ def test_model_solver_dae(self):
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
-    def test_model_solver_dae_bad_ics(self):
-        # Create model
+    def test_model_solver_dae_bad_ics_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var1 = pybamm.Variable("var1", domain=whole_cell)
         var2 = pybamm.Variable("var2", domain=whole_cell)
@@ -523,9 +523,9 @@ def test_model_solver_dae_bad_ics(self):
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
-    def test_model_solver_dae_events(self):
-        # Create model
+    def test_model_solver_dae_events_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var1 = pybamm.Variable("var1", domain=whole_cell)
         var2 = pybamm.Variable("var2", domain=whole_cell)
@@ -548,9 +548,9 @@ def test_model_solver_dae_events(self):
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
-    def test_model_solver_dae_with_jacobian(self):
-        # Create simple test model
+    def test_model_solver_dae_with_jacobian_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var1 = pybamm.Variable("var1", domain=whole_cell)
         var2 = pybamm.Variable("var2", domain=whole_cell)
@@ -584,9 +584,9 @@ def jacobian(t, y):
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
-    def test_solve_ode_model_with_dae_solver(self):
-        # Create model
+    def test_solve_ode_model_with_dae_solver_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         var = pybamm.Variable("var")
         model.rhs = {var: 0.1 * var}
         model.initial_conditions = {var: 1}
@@ -600,9 +600,9 @@ def test_solve_ode_model_with_dae_solver(self):
         np.testing.assert_array_equal(solution.t, t_eval)
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
 
-    def test_model_step_ode(self):
-        # Create model
+    def test_model_step_ode_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=whole_cell)
         model.rhs = {var: 0.1 * var}
@@ -629,9 +629,9 @@ def test_model_step_ode(self):
         concatenated_steps = np.concatenate((step_sol.y[0], step_sol_2.y[0, 1:]))
         np.testing.assert_allclose(solution.y[0], concatenated_steps)
 
-    def test_model_step_dae(self):
-        # Create model
+    def test_model_step_dae_python(self):
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var1 = pybamm.Variable("var1", domain=whole_cell)
         var2 = pybamm.Variable("var2", domain=whole_cell)
@@ -717,6 +717,22 @@ def test_model_solver_dae_events_casadi(self):
             np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
             np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
+    def test_solve_ode_model_with_dae_solver_casadi(self):
+        model = pybamm.BaseModel()
+        model.convert_to_format = "casadi"
+        var = pybamm.Variable("var")
+        model.rhs = {var: 0.1 * var}
+        model.initial_conditions = {var: 1}
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        t_eval = np.linspace(0, 1, 100)
+        solution = solver.solve(model, t_eval)
+        np.testing.assert_array_equal(solution.t, t_eval)
+        np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

tests/unit/test_solvers/test_scipy_solver.py

@@ -132,9 +132,10 @@ def jacobian(t, y):
             np.exp(1 + solution.t - np.exp(solution.t)), solution.y[1], rtol=1e-4
         )
 
-    def test_model_solver(self):
+    def test_model_solver_python(self):
         # Create model
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)
         model.rhs = {var: 0.1 * var}
@@ -158,9 +159,10 @@ def test_model_solver(self):
             solution.total_time, solution.solve_time + solution.set_up_time
         )
 
-    def test_model_solver_with_event(self):
+    def test_model_solver_with_event_python(self):
         # Create model
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)
         model.rhs = {var: -0.1 * var}
@@ -181,9 +183,10 @@ def test_model_solver_with_event(self):
         np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
-    def test_model_solver_ode_with_jacobian(self):
+    def test_model_solver_ode_with_jacobian_python(self):
         # Create model
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var1 = pybamm.Variable("var1", domain=whole_cell)
         var2 = pybamm.Variable("var2", domain=whole_cell)
@@ -232,9 +235,10 @@ def jacobian(t, y):
             np.ones((N, T.size)) * (T[np.newaxis, :] - np.exp(T[np.newaxis, :])),
         )
 
-    def test_model_step(self):
+    def test_model_step_python(self):
         # Create model
         model = pybamm.BaseModel()
+        model.convert_to_format = "python"
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)
         model.rhs = {var: 0.1 * var}