diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 98dc66eb3c..b45e1393d1 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -237,7 +237,7 @@ def version(formatted=False): from .solvers.scikits_ode_solver import ScikitsOdeSolver from .solvers.scikits_ode_solver import have_scikits_odes from .solvers.algebraic_solver import AlgebraicSolver -from .solvers.klu_sparse_solver import KLU +from .solvers.klu_sparse_solver import KLU, have_klu # diff --git a/pybamm/solvers/klu_sparse_solver.py b/pybamm/solvers/klu_sparse_solver.py index bdd72e9eb6..019d5ddff7 100644 --- a/pybamm/solvers/klu_sparse_solver.py +++ b/pybamm/solvers/klu_sparse_solver.py @@ -2,12 +2,20 @@ # Solver class using sundials with the KLU sparse linear solver # import pybamm - import numpy as np - -import klu import scipy.sparse as sparse +import importlib + +klu_spec = importlib.util.find_spec("klu") +if klu_spec is not None: + klu = importlib.util.module_from_spec(klu_spec) + klu_spec.loader.exec_module(klu) + + +def have_klu(): + return klu_spec is None + class KLU(pybamm.DaeSolver): """Solve a discretised model, using sundials with the KLU sparse linear solver. @@ -31,6 +39,8 @@ class KLU(pybamm.DaeSolver): def __init__( self, method="ida", tol=1e-8, root_method="lm", root_tol=1e-6, max_steps=1000 ): + if klu_spec is None: + raise ImportError("KLU is not installed") super().__init__(method, tol, root_method, root_tol, max_steps) diff --git a/pybamm/solvers/scikits_ode_solver.py b/pybamm/solvers/scikits_ode_solver.py index 17cfff62a0..bbeed30ccd 100644 --- a/pybamm/solvers/scikits_ode_solver.py +++ b/pybamm/solvers/scikits_ode_solver.py @@ -134,7 +134,7 @@ def jac_times_setupfn(t, y, fy, userdata): np.transpose(sol.values.y), sol.roots.t, np.transpose(sol.roots.y), - termination + termination, ) else: raise pybamm.SolverError(sol.message) diff --git a/tests/unit/test_solvers/test_klu_sparse_solver.py b/tests/unit/test_solvers/test_klu_sparse_solver.py index 4db3b2622c..c48f7340f1 100644 --- a/tests/unit/test_solvers/test_klu_sparse_solver.py +++ b/tests/unit/test_solvers/test_klu_sparse_solver.py @@ -8,7 +8,7 @@ from tests import get_mesh_for_testing, get_discretisation_for_testing -@unittest.skipIf(pybamm.have_scikits_odes(), "scikits.odes not installed") +@unittest.skipIf(pybamm.have_klu(), "klu is not installed") class TestKLUSolver(unittest.TestCase): def test_ida_roberts_klu(self): # this test implements a python version of the ida Roberts @@ -28,12 +28,6 @@ def jac(t, y): J[1][1] = -1.0 return sparse.csr_matrix(J) - # def events(t, y): - # g = np.zeros((num_of_events,)) - # g[0] = y[0] - 1.0 - # g[1] = y[1] - 0.0 - # return g - def event_1(t, y): return y[0] - 0.2 @@ -67,596 +61,21 @@ def alg(t, y): solution = solver.integrate(res, y0, t_eval, events, mass_matrix, jac) - print("hello") - - # test that final time is triggered by event + # test that final time is time of event + # y = 0.1 t + y0 so y=0.2 when t=2 + np.testing.assert_array_almost_equal(solution.t[-1], 2.0) # test that final value is the event value + np.testing.assert_array_almost_equal(solution.y[0, -1], 0.2) - # test that y[2] remains constant - - # test that y[1] = to true solution - - # def test_ode_integrate_with_jacobian(self): - # # Linear - # solver = pybamm.KLU(tol=1e-8) - - # def linear_ode(t, y): - # return np.array([0.5, 2 - y[0]]) - - # J = np.array([[0.0, 0.0], [-1.0, 0.0]]) - # sJ = sparse.csr_matrix(J) - - # def jacobian(t, y): - # return J - - # def sparse_jacobian(t, y): - # return sJ - - # y0 = np.array([0.0, 0.0]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate(linear_ode, y0, t_eval, jacobian=jacobian) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - # np.testing.assert_allclose( - # 2.0 * solution.t - 0.25 * solution.t ** 2, solution.y[1], rtol=1e-4 - # ) - - # y0 = np.array([0.0, 0.0]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate(linear_ode, y0, t_eval, jacobian=sparse_jacobian) - - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose( - # 2.0 * solution.t - 0.25 * solution.t ** 2, solution.y[1], rtol=1e-4 - # ) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - - # solver = pybamm.KLU(tol=1e-8) - - # solution = solver.integrate(linear_ode, y0, t_eval, jacobian=jacobian) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - # np.testing.assert_allclose( - # 2.0 * solution.t - 0.25 * solution.t ** 2, solution.y[1], rtol=1e-4 - # ) - - # solution = solver.integrate(linear_ode, y0, t_eval, jacobian=sparse_jacobian) - - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose( - # 2.0 * solution.t - 0.25 * solution.t ** 2, solution.y[1], rtol=1e-4 - # ) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - - # # Nonlinear exponential grwoth - # solver = pybamm.KLU(tol=1e-8) - - # def exponential_growth(t, y): - # return np.array([y[0], (1.0 - y[0]) * y[1]]) - - # def jacobian(t, y): - # return np.array([[1.0, 0.0], [-y[1], 1 - y[0]]]) - - # def sparse_jacobian(t, y): - # return sparse.csr_matrix(jacobian(t, y)) - - # y0 = np.array([1.0, 1.0]) - # t_eval = np.linspace(0, 1, 100) - - # solution = solver.integrate(exponential_growth, y0, t_eval, jacobian=jacobian) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(np.exp(solution.t), solution.y[0], rtol=1e-4) - # np.testing.assert_allclose( - # np.exp(1 + solution.t - np.exp(solution.t)), solution.y[1], rtol=1e-4 - # ) - - # solution = solver.integrate( - # exponential_growth, y0, t_eval, jacobian=sparse_jacobian - # ) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(np.exp(solution.t), solution.y[0], rtol=1e-4) - # np.testing.assert_allclose( - # np.exp(1 + solution.t - np.exp(solution.t)), solution.y[1], rtol=1e-4 - # ) - - # solver = pybamm.KLU(tol=1e-8, linsolver="spgmr") - - # solution = solver.integrate(exponential_growth, y0, t_eval, jacobian=jacobian) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(np.exp(solution.t), solution.y[0], rtol=1e-4) - # np.testing.assert_allclose( - # np.exp(1 + solution.t - np.exp(solution.t)), solution.y[1], rtol=1e-4 - # ) - - # solution = solver.integrate( - # exponential_growth, y0, t_eval, jacobian=sparse_jacobian - # ) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(np.exp(solution.t), solution.y[0], rtol=1e-4) - # np.testing.assert_allclose( - # np.exp(1 + solution.t - np.exp(solution.t)), solution.y[1], rtol=1e-4 - # ) - - # def test_dae_integrate(self): - # # Constant - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def constant_growth_dae(t, y, ydot): - # return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]] - - # y0 = np.array([0, 0]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate(constant_growth_dae, y0, t_eval) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - # np.testing.assert_allclose(1.0 * solution.t, solution.y[1]) - - # # Exponential decay - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def exponential_decay_dae(t, y, ydot): - # return [-0.1 * y[0] - ydot[0], 2 * y[0] - y[1]] - - # y0 = np.array([1, 2]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate(exponential_decay_dae, y0, t_eval) - # 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)) - # self.assertEqual(solution.termination, "final time") - - # def test_dae_integrate_failure(self): - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def constant_growth_dae(t, y, ydot): - # return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]] - - # y0 = np.array([0, 1]) - # t_eval = np.linspace(0, 1, 100) - # with self.assertRaises(pybamm.SolverError): - # solver.integrate(constant_growth_dae, y0, t_eval) - - # def test_dae_integrate_bad_ics(self): - # # Constant - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def constant_growth_dae(t, y, ydot): - # return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]] - - # def constant_growth_dae_rhs(t, y): - # return np.array([constant_growth_dae(t, y, [0])[0]]) - - # def constant_growth_dae_algebraic(t, y): - # return np.array([constant_growth_dae(t, y, [0])[1]]) - - # y0_guess = np.array([0, 1]) - # t_eval = np.linspace(0, 1, 100) - # y0 = solver.calculate_consistent_initial_conditions( - # constant_growth_dae_rhs, constant_growth_dae_algebraic, y0_guess - # ) - # # check y0 - # np.testing.assert_array_equal(y0, [0, 0]) - # # check dae solutions - # solution = solver.integrate(constant_growth_dae, y0, t_eval) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - # np.testing.assert_allclose(1.0 * solution.t, solution.y[1]) - - # # Exponential decay - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def exponential_decay_dae(t, y, ydot): - # return [-0.1 * y[0] - ydot[0], 2 * y[0] - y[1]] - - # y0 = np.array([1, 2]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate(exponential_decay_dae, y0, t_eval) - # 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_dae_integrate_with_event(self): - # # Constant - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def constant_growth_dae(t, y, ydot): - # return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]] - - # def y0_eq_2(t, y): - # return y[0] - 2 - - # def y1_eq_5(t, y): - # return y[1] - 5 - - # y0 = np.array([0, 0]) - # t_eval = np.linspace(0, 7, 100) - # solution = solver.integrate( - # constant_growth_dae, y0, t_eval, events=[y0_eq_2, y1_eq_5] - # ) - # self.assertLess(len(solution.t), len(t_eval)) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - # np.testing.assert_allclose(1.0 * solution.t, solution.y[1]) - # np.testing.assert_array_less(solution.y[0], 2) - # np.testing.assert_array_less(solution.y[1], 5) - # np.testing.assert_allclose(solution.t_event, 4) - # np.testing.assert_allclose(solution.y_event[0], 2) - # np.testing.assert_allclose(solution.y_event[1], 4) - # self.assertEqual(solution.termination, "event") - - # # Exponential decay - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def exponential_decay_dae(t, y, ydot): - # return np.array([-0.1 * y[0] - ydot[0], 2 * y[0] - y[1]]) - - # def y0_eq_0pt9(t, y): - # return y[0] - 0.9 - - # def t_eq_0pt5(t, y): - # return t - 0.5 - - # y0 = np.array([1, 2]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate( - # exponential_decay_dae, y0, t_eval, events=[y0_eq_0pt9, t_eq_0pt5] - # ) - - # self.assertLess(len(solution.t), len(t_eval)) - # 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)) - # np.testing.assert_array_less(0.9, solution.y[0]) - # np.testing.assert_array_less(solution.t, 0.5) - # np.testing.assert_allclose(solution.t_event, 0.5) - # np.testing.assert_allclose(solution.y_event[0], np.exp(-0.1 * 0.5)) - # np.testing.assert_allclose(solution.y_event[1], 2 * np.exp(-0.1 * 0.5)) - # self.assertEqual(solution.termination, "event") - - # def test_dae_integrate_with_jacobian(self): - # # Constant - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def constant_growth_dae(t, y, ydot): - # return np.array([0.5 * np.ones_like(y[0]) - ydot[0], 2.0 * y[0] - y[1]]) - - # mass_matrix = np.array([[1.0, 0.0], [0.0, 0.0]]) - - # def jacobian(t, y): - # return np.array([[0.0, 0.0], [2.0, -1.0]]) - - # y0 = np.array([0.0, 0.0]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate( - # constant_growth_dae, y0, t_eval, mass_matrix=mass_matrix, jacobian=jacobian - # ) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - # np.testing.assert_allclose(1.0 * solution.t, solution.y[1]) - - # # Nonlinear (tests when Jacobian a function of y) - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def nonlinear_dae(t, y, ydot): - # return np.array([0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] ** 2 - y[1]]) - - # mass_matrix = np.array([[1.0, 0.0], [0.0, 0.0]]) - - # def jacobian(t, y): - # return np.array([[0.0, 0.0], [4.0 * y[0], -1.0]]) - - # y0 = np.array([0.0, 0.0]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate( - # nonlinear_dae, y0, t_eval, mass_matrix=mass_matrix, jacobian=jacobian - # ) - # np.testing.assert_array_equal(solution.t, t_eval) - # np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) - # np.testing.assert_allclose(0.5 * solution.t ** 2, solution.y[1]) - - # def test_dae_integrate_with_non_unity_mass(self): - # # Constant - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # def constant_growth_dae(t, y, ydot): - # return np.array([0.5 * np.ones_like(y[0]) - 4 * ydot[0], 2.0 * y[0] - y[1]]) - - # mass_matrix = np.array([[4.0, 0.0], [0.0, 0.0]]) - - # def jacobian(t, y): - # return np.array([[0.0, 0.0], [2.0, -1.0]]) - - # y0 = np.array([0.0, 0.0]) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.integrate( - # constant_growth_dae, y0, t_eval, mass_matrix=mass_matrix, jacobian=jacobian - # ) - # np.testing.assert_array_equal(solution.t, t_eval) - # 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 - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var = pybamm.Variable("var", domain=whole_cell) - # model.rhs = {var: 0.1 * var} - # model.initial_conditions = {var: 1} - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # # Solve - # solver = pybamm.KLU(tol=1e-9) - # 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)) - - # # Test time - # self.assertGreater( - # solution.total_time, solution.solve_time + solution.set_up_time - # ) - - # def test_model_solver_ode_events(self): - # # Create model - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var = pybamm.Variable("var", domain=whole_cell) - # model.rhs = {var: 0.1 * var} - # model.initial_conditions = {var: 1} - # model.events = { - # "2 * var = 2.5": pybamm.min(2 * var - 2.5), - # "var = 1.5": pybamm.min(var - 1.5), - # } - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # # Solve - # solver = pybamm.KLU(tol=1e-9) - # t_eval = np.linspace(0, 10, 100) - # solution = solver.solve(model, t_eval) - # np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t)) - # 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 - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var1 = pybamm.Variable("var1", domain=whole_cell) - # var2 = pybamm.Variable("var2", domain=whole_cell) - # model.rhs = {var1: var1, var2: 1 - var1} - # model.initial_conditions = {var1: 1.0, var2: -1.0} - # model.variables = {"var1": var1, "var2": var2} - - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # # Add user-supplied Jacobian to model - # mesh = get_mesh_for_testing() - # combined_submesh = mesh.combine_submeshes( - # "negative electrode", "separator", "positive electrode" - # ) - # N = combined_submesh[0].npts - - # # construct jacobian in order of model.rhs - # J = [] - # for var in model.rhs.keys(): - # if var.id == var1.id: - # J.append([np.eye(N), np.zeros((N, N))]) - # else: - # J.append([-1.0 * np.eye(N), np.zeros((N, N))]) - - # J = np.block(J) - - # def jacobian(t, y): - # return J - - # model.jacobian = jacobian - - # # Solve - # solver = pybamm.KLU(tol=1e-9) - # t_eval = np.linspace(0, 1, 100) - # solution = solver.solve(model, t_eval) - # np.testing.assert_array_equal(solution.t, t_eval) - - # T, Y = solution.t, solution.y - # np.testing.assert_array_almost_equal( - # model.variables["var1"].evaluate(T, Y), - # np.ones((N, T.size)) * np.exp(T[np.newaxis, :]), - # ) - # np.testing.assert_array_almost_equal( - # model.variables["var2"].evaluate(T, Y), - # np.ones((N, T.size)) * (T[np.newaxis, :] - np.exp(T[np.newaxis, :])), - # ) - - # def test_model_solver_dae(self): - # # Create model - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var1 = pybamm.Variable("var1", domain=whole_cell) - # var2 = pybamm.Variable("var2", domain=whole_cell) - # model.rhs = {var1: 0.1 * var1} - # model.algebraic = {var2: 2 * var1 - var2} - # model.initial_conditions = {var1: 1, var2: 2} - # model.use_jacobian = False - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # # Solve - # solver = pybamm.ScikitsDaeSolver(tol=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)) - # np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t)) - - # # Test time - # self.assertGreater( - # solution.total_time, solution.solve_time + solution.set_up_time - # ) - - # def test_model_solver_dae_bad_ics(self): - # # Create model - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var1 = pybamm.Variable("var1", domain=whole_cell) - # var2 = pybamm.Variable("var2", domain=whole_cell) - # model.rhs = {var1: 0.1 * var1} - # model.algebraic = {var2: 2 * var1 - var2} - # model.initial_conditions = {var1: 1, var2: 3} - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # # Solve - # solver = pybamm.ScikitsDaeSolver(tol=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)) - # np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t)) - - # def test_model_solver_dae_events(self): - # # Create model - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var1 = pybamm.Variable("var1", domain=whole_cell) - # var2 = pybamm.Variable("var2", domain=whole_cell) - # model.rhs = {var1: 0.1 * var1} - # model.algebraic = {var2: 2 * var1 - var2} - # model.initial_conditions = {var1: 1, var2: 2} - # model.events = { - # "var1 = 1.5": pybamm.min(var1 - 1.5), - # "var2 = 2.5": pybamm.min(var2 - 2.5), - # } - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # # Solve - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - # t_eval = np.linspace(0, 5, 100) - # solution = solver.solve(model, t_eval) - # np.testing.assert_array_less(solution.y[0], 1.5) - # np.testing.assert_array_less(solution.y[-1], 2.5) - # 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 - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var1 = pybamm.Variable("var1", domain=whole_cell) - # var2 = pybamm.Variable("var2", domain=whole_cell) - # model.rhs = {var1: 0.1 * var1} - # model.algebraic = {var2: 2 * var1 - var2} - # model.initial_conditions = {var1: 1.0, var2: 2.0} - # model.initial_conditions_ydot = {var1: 0.1, var2: 0.2} - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # # Add user-supplied Jacobian to model - # mesh = get_mesh_for_testing() - # combined_submesh = mesh.combine_submeshes( - # "negative electrode", "separator", "positive electrode" - # ) - # N = combined_submesh[0].npts - - # def jacobian(t, y): - # return np.block( - # [ - # [0.1 * np.eye(N), np.zeros((N, N))], - # [2.0 * np.eye(N), -1.0 * np.eye(N)], - # ] - # ) - - # model.jacobian = jacobian - - # # Solve - # solver = pybamm.ScikitsDaeSolver(tol=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)) - # 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 - # model = pybamm.BaseModel() - # 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(tol=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)) - - # def test_model_step_ode(self): - # # Create model - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var = pybamm.Variable("var", domain=whole_cell) - # model.rhs = {var: 0.1 * var} - # model.initial_conditions = {var: 1} - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # solver = pybamm.KLU(tol=1e-9) - - # # Step once - # dt = 0.1 - # step_sol = solver.step(model, dt) - # np.testing.assert_array_equal(step_sol.t, [0, dt]) - # np.testing.assert_allclose(step_sol.y[0], np.exp(0.1 * step_sol.t)) - - # # Step again (return 5 points) - # step_sol_2 = solver.step(model, dt, npts=5) - # np.testing.assert_array_equal(step_sol_2.t, np.linspace(dt, 2 * dt, 5)) - # np.testing.assert_allclose(step_sol_2.y[0], np.exp(0.1 * step_sol_2.t)) - - # # Check steps give same solution as solve - # t_eval = np.concatenate((step_sol.t, step_sol_2.t[1:])) - # solution = solver.solve(model, t_eval) - # 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 - # model = pybamm.BaseModel() - # whole_cell = ["negative electrode", "separator", "positive electrode"] - # var1 = pybamm.Variable("var1", domain=whole_cell) - # var2 = pybamm.Variable("var2", domain=whole_cell) - # model.rhs = {var1: 0.1 * var1} - # model.algebraic = {var2: 2 * var1 - var2} - # model.initial_conditions = {var1: 1, var2: 2} - # model.use_jacobian = False - # disc = get_discretisation_for_testing() - # disc.process_model(model) - - # solver = pybamm.ScikitsDaeSolver(tol=1e-8) - - # # Step once - # dt = 0.1 - # step_sol = solver.step(model, dt) - # np.testing.assert_array_equal(step_sol.t, [0, dt]) - # np.testing.assert_allclose(step_sol.y[0], np.exp(0.1 * step_sol.t)) - # np.testing.assert_allclose(step_sol.y[-1], 2 * np.exp(0.1 * step_sol.t)) - - # # Step again (return 5 points) - # step_sol_2 = solver.step(model, dt, npts=5) - # np.testing.assert_array_equal(step_sol_2.t, np.linspace(dt, 2 * dt, 5)) - # np.testing.assert_allclose(step_sol_2.y[0], np.exp(0.1 * step_sol_2.t)) - # np.testing.assert_allclose(step_sol_2.y[-1], 2 * np.exp(0.1 * step_sol_2.t)) - - # # append solutions - # step_sol.append(step_sol_2) + # test that y[1] remains constant + np.testing.assert_array_almost_equal( + solution.y[1, :], np.ones(solution.t.shape) + ) - # # Check steps give same solution as solve - # t_eval = step_sol.t - # solution = solver.solve(model, t_eval) - # np.testing.assert_allclose(solution.y[0], step_sol.y[0, :]) - # np.testing.assert_allclose(solution.y[-1], step_sol.y[-1, :]) + # test that y[0] = to true solution + true_solution = 0.1 * solution.t + np.testing.assert_array_almost_equal(solution.y[0, :], true_solution) if __name__ == "__main__":