diff --git a/CHANGELOG.md b/CHANGELOG.md index 8c2c65ce0b..b2e64761ed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) +## Bug fixes + +- For simulations with events that cause the simulation to stop early, the sensitivities could be evaluated incorrectly to zero ([#2331](https://github.com/pybamm-team/PyBaMM/pull/2337)) + # [v22.9](https://github.com/pybamm-team/PyBaMM/tree/v22.9) - 2022-09-30 ## Features diff --git a/pybamm/solvers/c_solvers/idaklu.cpp b/pybamm/solvers/c_solvers/idaklu.cpp index e23f037f9a..9241487419 100644 --- a/pybamm/solvers/c_solvers/idaklu.cpp +++ b/pybamm/solvers/c_solvers/idaklu.cpp @@ -283,7 +283,8 @@ Solution solve_python(np_array t_np, np_array y0_np, np_array yp0_np, void *ida_mem; // pointer to memory N_Vector yy, yp, avtol; // y, y', and absolute tolerance N_Vector *yyS, *ypS; // y, y' for sensitivities - realtype rtol, *yval, *ypval, *atval, *ySval; + realtype rtol, *yval, *ypval, *atval; + std::vector ySval(number_of_parameters); int retval; SUNMatrix J; SUNLinearSolver LS; @@ -300,9 +301,6 @@ Solution solve_python(np_array t_np, np_array y0_np, np_array yp0_np, // set initial value yval = N_VGetArrayPointer(yy); - if (number_of_parameters > 0) { - ySval = N_VGetArrayPointer(yyS[0]); - } ypval = N_VGetArrayPointer(yp); atval = N_VGetArrayPointer(avtol); int i; @@ -314,6 +312,7 @@ Solution solve_python(np_array t_np, np_array y0_np, np_array yp0_np, } for (int is = 0 ; is < number_of_parameters; is++) { + ySval[is] = N_VGetArrayPointer(yyS[is]); N_VConst(RCONST(0.0), yyS[is]); N_VConst(RCONST(0.0), ypS[is]); } @@ -376,7 +375,7 @@ Solution solve_python(np_array t_np, np_array y0_np, np_array yp0_np, for (int j = 0; j < number_of_parameters; j++) { const int base_index = j * number_of_timesteps * number_of_states; for (int k = 0; k < number_of_states; k++) { - yS_return[base_index + k] = ySval[j * number_of_states + k]; + yS_return[base_index + k] = ySval[j][k]; } } @@ -417,7 +416,7 @@ Solution solve_python(np_array t_np, np_array y0_np, np_array yp0_np, const int base_index = j * number_of_timesteps * number_of_states + t_i * number_of_states; for (int k = 0; k < number_of_states; k++) { - yS_return[base_index + k] = ySval[j * number_of_states + k]; + yS_return[base_index + k] = ySval[j][k]; } } t_i += 1; @@ -445,7 +444,7 @@ Solution solve_python(np_array t_np, np_array y0_np, np_array yp0_np, np_array t_ret = np_array(t_i, &t_return[0]); np_array y_ret = np_array(t_i * number_of_states, &y_return[0]); np_array yS_ret = np_array( - std::vector{number_of_parameters, t_i, number_of_states}, + std::vector{number_of_parameters, number_of_timesteps, number_of_states}, &yS_return[0] ); diff --git a/pybamm/solvers/c_solvers/idaklu_casadi.cpp b/pybamm/solvers/c_solvers/idaklu_casadi.cpp index 0f123df4a0..f611c6ebc3 100644 --- a/pybamm/solvers/c_solvers/idaklu_casadi.cpp +++ b/pybamm/solvers/c_solvers/idaklu_casadi.cpp @@ -15,15 +15,6 @@ class CasadiFunction { size_t sz_iw; size_t sz_w; m_func.sz_work(sz_arg, sz_res, sz_iw, sz_w); - //std::cout << "name = "<< m_func.name() << " arg = " << sz_arg << " res = " << sz_res << " iw = " << sz_iw << " w = " << sz_w << std::endl; - //for (int i = 0; i < sz_arg; i++) { - // std::cout << "Sparsity for input " << i << std::endl; - // const Sparsity& sparsity = m_func.sparsity_in(i); - //} - //for (int i = 0; i < sz_res; i++) { - // std::cout << "Sparsity for output " << i << std::endl; - // const Sparsity& sparsity = m_func.sparsity_out(i); - //} m_arg.resize(sz_arg); m_res.resize(sz_res); m_iw.resize(sz_iw); @@ -104,17 +95,6 @@ int residual_casadi(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr, CasadiFunctions *p_python_functions = static_cast(user_data); - //std::cout << "RESIDUAL t = " << tres << " y = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(yy)[i] << " "; - //} - //std::cout << "] yp = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(yp)[i] << " "; - //} - //std::cout << "]" << std::endl; - // args are t, y, put result in rr - py::buffer_info input_buf = p_python_functions->inputs.request(); p_python_functions->rhs_alg.m_arg[0] = &tres; p_python_functions->rhs_alg.m_arg[1] = NV_DATA_S(yy); @@ -122,41 +102,18 @@ int residual_casadi(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr, p_python_functions->rhs_alg.m_res[0] = NV_DATA_S(rr); p_python_functions->rhs_alg(); - //std::cout << "rhs_alg = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(rr)[i] << " "; - //} - //std::cout << "]" << std::endl; - realtype *tmp = p_python_functions->get_tmp(); - //std::cout << "tmp before = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << tmp[i] << " "; - //} - //std::cout << "]" << std::endl; + // args is yp, put result in tmp p_python_functions->mass_action.m_arg[0] = NV_DATA_S(yp); p_python_functions->mass_action.m_res[0] = tmp; p_python_functions->mass_action(); - //std::cout << "tmp = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << tmp[i] << " "; - //} - //std::cout << "]" << std::endl; - // AXPY: y <- a*x + y const int ns = p_python_functions->number_of_states; casadi_axpy(ns, -1., tmp, NV_DATA_S(rr)); - //std::cout << "residual = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(rr)[i] << " "; - //} - //std::cout << "]" << std::endl; - // now rr has rhs_alg(t, y) - mass_matrix * yp - return 0; } @@ -246,16 +203,9 @@ int jacobian_casadi(realtype tt, realtype cj, N_Vector yy, N_Vector yp, const int n_row_vals = jac_times_cjmass_rowvals.request().size; auto p_jac_times_cjmass_rowvals = jac_times_cjmass_rowvals.unchecked<1>(); - //std::cout << "jac_data = ["; - //for (int i = 0; i < p_python_functions->number_of_nnz; i++) { - // std::cout << jac_data[i] << " "; - //} - //std::cout << "]" << std::endl; - // just copy across row vals (do I need to do this every time?) // (or just in the setup?) for (int i = 0; i < n_row_vals; i++) { - //std::cout << "check row vals " << jac_rowvals[i] << " " << p_jac_times_cjmass_rowvals[i] << std::endl; jac_rowvals[i] = p_jac_times_cjmass_rowvals[i]; } @@ -265,7 +215,6 @@ int jacobian_casadi(realtype tt, realtype cj, N_Vector yy, N_Vector yp, // just copy across col ptrs (do I need to do this every time?) for (int i = 0; i < n_col_ptrs; i++) { - //std::cout << "check col ptrs " << jac_colptrs[i] << " " << p_jac_times_cjmass_colptrs[i] << std::endl; jac_colptrs[i] = p_jac_times_cjmass_colptrs[i]; } @@ -278,17 +227,6 @@ int events_casadi(realtype t, N_Vector yy, N_Vector yp, realtype *events_ptr, CasadiFunctions *p_python_functions = static_cast(user_data); - //std::cout << "EVENTS" << std::endl; - //std::cout << "t = " << t << " y = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(yy)[i] << " "; - //} - //std::cout << "] yp = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(yp)[i] << " "; - //} - //std::cout << "]" << std::endl; - // args are t, y, put result in events_ptr py::buffer_info input_buf = p_python_functions->inputs.request(); p_python_functions->events.m_arg[0] = &t; @@ -297,12 +235,6 @@ int events_casadi(realtype t, N_Vector yy, N_Vector yp, realtype *events_ptr, p_python_functions->events.m_res[0] = events_ptr; p_python_functions->events(); - //std::cout << "events = ["; - //for (int i = 0; i < p_python_functions->number_of_events; i++) { - // std::cout << events_ptr[i] << " "; - //} - //std::cout << "]" << std::endl; - return (0); } @@ -336,32 +268,6 @@ int sensitivities_casadi(int Ns, realtype t, N_Vector yy, N_Vector yp, const int np = p_python_functions->number_of_parameters; - //std::cout << "SENS t = " << t << " y = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(yy)[i] << " "; - //} - //std::cout << "] yp = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(yp)[i] << " "; - //} - //std::cout << "] yS = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(yS[0])[i] << " "; - //} - //std::cout << "] ypS = ["; - //for (int i = 0; i < p_python_functions->number_of_states; i++) { - // std::cout << NV_DATA_S(ypS[0])[i] << " "; - //} - //std::cout << "]" << std::endl; - - //for (int i = 0; i < np; i++) { - // std::cout << "dF/dp before = [" << i << "] = ["; - // for (int j = 0; j < p_python_functions->number_of_states; j++) { - // std::cout << NV_DATA_S(resvalS[i])[j] << " "; - // } - // std::cout << "]" << std::endl; - //} - // args are t, y put result in rr py::buffer_info input_buf = p_python_functions->inputs.request(); p_python_functions->sens.m_arg[0] = &t; @@ -374,13 +280,6 @@ int sensitivities_casadi(int Ns, realtype t, N_Vector yy, N_Vector yp, p_python_functions->sens(); for (int i = 0; i < np; i++) { - //std::cout << "dF/dp = [" << i << "] = ["; - //for (int j = 0; j < p_python_functions->number_of_states; j++) { - // std::cout << NV_DATA_S(resvalS[i])[j] << " "; - //} - //std::cout << "]" << std::endl; - - // put (∂F/∂y)s i (t) in tmp realtype *tmp = p_python_functions->get_tmp(); p_python_functions->jac_action.m_arg[0] = &t; @@ -390,12 +289,6 @@ int sensitivities_casadi(int Ns, realtype t, N_Vector yy, N_Vector yp, p_python_functions->jac_action.m_res[0] = tmp; p_python_functions->jac_action(); - //std::cout << "jac_action = [" << i << "] = ["; - //for (int j = 0; j < p_python_functions->number_of_states; j++) { - // std::cout << tmp[j] << " "; - //} - //std::cout << "]" << std::endl; - const int ns = p_python_functions->number_of_states; casadi_axpy(ns, 1., tmp, NV_DATA_S(resvalS[i])); @@ -404,22 +297,9 @@ int sensitivities_casadi(int Ns, realtype t, N_Vector yy, N_Vector yp, p_python_functions->mass_action.m_res[0] = tmp; p_python_functions->mass_action(); - //std::cout << "mass_Action = [" << i << "] = ["; - //for (int j = 0; j < p_python_functions->number_of_states; j++) { - // std::cout << tmp[j] << " "; - //} - //std::cout << "]" << std::endl; - - // (∂F/∂y)s i (t)+(∂F/∂ ẏ) ṡ i (t)+(∂F/∂p i ) // AXPY: y <- a*x + y casadi_axpy(ns, -1., tmp, NV_DATA_S(resvalS[i])); - - //std::cout << "resvalS[" << i << "] = ["; - //for (int j = 0; j < p_python_functions->number_of_states; j++) { - // std::cout << NV_DATA_S(resvalS[i])[j] << " "; - //} - //std::cout << "]" << std::endl; } return 0; @@ -457,7 +337,8 @@ Solution solve_casadi(np_array t_np, np_array y0_np, np_array yp0_np, void *ida_mem; // pointer to memory N_Vector yy, yp, avtol; // y, y', and absolute tolerance N_Vector *yyS, *ypS; // y, y' for sensitivities - realtype rtol, *yval, *ypval, *atval, *ySval; + realtype rtol, *yval, *ypval, *atval; + std::vector ySval(number_of_parameters); int retval; SUNMatrix J; SUNLinearSolver LS; @@ -474,9 +355,6 @@ Solution solve_casadi(np_array t_np, np_array y0_np, np_array yp0_np, // set initial value yval = N_VGetArrayPointer(yy); - if (number_of_parameters > 0) { - ySval = N_VGetArrayPointer(yyS[0]); - } ypval = N_VGetArrayPointer(yp); atval = N_VGetArrayPointer(avtol); int i; @@ -488,6 +366,7 @@ Solution solve_casadi(np_array t_np, np_array y0_np, np_array yp0_np, } for (int is = 0 ; is < number_of_parameters; is++) { + ySval[is] = N_VGetArrayPointer(yyS[is]); N_VConst(RCONST(0.0), yyS[is]); N_VConst(RCONST(0.0), ypS[is]); } @@ -547,37 +426,6 @@ Solution solve_casadi(np_array t_np, np_array y0_np, np_array yp0_np, SUNLinSolInitialize(LS); - //std::cout << "number of jacobian nz = " << jac_times_cjmass_nnz << std::endl; - //N_Vector rr, tmp1, tmp2, tmp3; - //rr = N_VNew_Serial(number_of_states); - //tmp1 = N_VNew_Serial(number_of_states); - //tmp2 = N_VNew_Serial(number_of_states); - //tmp3 = N_VNew_Serial(number_of_states); - //std::vector events_vect(number_of_events); - //auto start = std::chrono::high_resolution_clock::now(); - //for (int i = 0; i < 10; i++) { - // residual_casadi(0, yy, yp, rr, user_data); - //} - //auto end = std::chrono::high_resolution_clock::now(); - //std::chrono::duration diff = end - start; - //std::cout << "Time to call residual " << diff.count() << " s\n"; - - //start = std::chrono::high_resolution_clock::now(); - //for (int i = 0; i < 10; i++) { - // events_casadi(0, yy, yp, events_vect.data(), user_data); - //} - //end = std::chrono::high_resolution_clock::now(); - //diff = end - start; - //std::cout << "Time to call events " << diff.count() << " s\n"; - - //start = std::chrono::high_resolution_clock::now(); - //for (int i = 0; i < 10; i++) { - // jacobian_casadi(0, 1, yy, yp, rr, J, user_data, tmp1, tmp2, tmp3); - //} - //end = std::chrono::high_resolution_clock::now(); - //diff = end - start; - //std::cout << "Time to call jacobian " << diff.count() << " s\n"; - int t_i = 1; realtype tret; realtype t_next; @@ -609,7 +457,7 @@ Solution solve_casadi(np_array t_np, np_array y0_np, np_array yp0_np, for (int j = 0; j < number_of_parameters; j++) { const int base_index = j * number_of_timesteps * number_of_states; for (int k = 0; k < number_of_states; k++) { - yS_return[base_index + k] = ySval[j * number_of_states + k]; + yS_return[base_index + k] = ySval[j][k]; } } @@ -649,7 +497,7 @@ Solution solve_casadi(np_array t_np, np_array y0_np, np_array yp0_np, const int base_index = j * number_of_timesteps * number_of_states + t_i * number_of_states; for (int k = 0; k < number_of_states; k++) { - yS_return[base_index + k] = ySval[j * number_of_states + k]; + yS_return[base_index + k] = ySval[j][k]; } } t_i += 1; @@ -665,7 +513,7 @@ Solution solve_casadi(np_array t_np, np_array y0_np, np_array yp0_np, np_array t_ret = np_array(t_i, &t_return[0], free_t_when_done); np_array y_ret = np_array(t_i * number_of_states, &y_return[0], free_y_when_done); np_array yS_ret = np_array( - std::vector{number_of_parameters, t_i, number_of_states}, + std::vector{number_of_parameters, number_of_timesteps, number_of_states}, &yS_return[0], free_yS_when_done ); @@ -729,10 +577,6 @@ Solution solve_casadi(np_array t_np, np_array y0_np, np_array yp0_np, IDAFree(&ida_mem); - //std::cout << "finished solving 9" << std::endl; - - //std::cout << "finished solving 10" << std::endl; - return sol; } diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index f823a20f60..ba4c810c58 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -246,6 +246,86 @@ def test_ida_roberts_klu_sensitivities(self): np.testing.assert_array_almost_equal(dyda_ida, dyda_fd) + def test_sensitivities_with_events(self): + # this test implements a python version of the ida Roberts + # example provided in sundials + # see sundials ida examples pdf + for form in ["python", "casadi", "jax"]: + if form == "jax" and not pybamm.have_jax(): + continue + if form == "casadi": + root_method = "casadi" + else: + root_method = "lm" + model = pybamm.BaseModel() + model.convert_to_format = form + u = pybamm.Variable("u") + v = pybamm.Variable("v") + a = pybamm.InputParameter("a") + b = pybamm.InputParameter("b") + model.rhs = {u: a * v + b} + model.algebraic = {v: 1 - v} + model.initial_conditions = {u: 0, v: 1} + model.events = [pybamm.Event("1", 0.2 - u)] + + disc = pybamm.Discretisation() + disc.process_model(model) + + solver = pybamm.IDAKLUSolver(root_method=root_method) + + t_eval = np.linspace(0, 3, 100) + a_value = 0.1 + b_value = 0.0 + + # solve first without sensitivities + sol = solver.solve( + model, + t_eval, + inputs={"a": a_value, "b": b_value}, + calculate_sensitivities=True, + ) + + # test that y[1] remains constant + np.testing.assert_array_almost_equal(sol.y[1, :], np.ones(sol.t.shape)) + + # test that y[0] = to true solution + true_solution = a_value * sol.t + np.testing.assert_array_almost_equal(sol.y[0, :], true_solution) + + # evaluate the sensitivities using idas + dyda_ida = sol.sensitivities["a"] + dydb_ida = sol.sensitivities["b"] + + # evaluate the sensitivities using finite difference + h = 1e-6 + sol_plus = solver.solve( + model, t_eval, inputs={"a": a_value + 0.5 * h, "b": b_value} + ) + sol_neg = solver.solve( + model, t_eval, inputs={"a": a_value - 0.5 * h, "b": b_value} + ) + max_index = min(sol_plus.y.shape[1], sol_neg.y.shape[1]) - 1 + dyda_fd = (sol_plus.y[:, :max_index] - sol_neg.y[:, :max_index]) / h + dyda_fd = dyda_fd.transpose().reshape(-1, 1) + + np.testing.assert_array_almost_equal( + dyda_ida[: (2 * max_index), :], dyda_fd + ) + + sol_plus = solver.solve( + model, t_eval, inputs={"a": a_value, "b": b_value + 0.5 * h} + ) + sol_neg = solver.solve( + model, t_eval, inputs={"a": a_value, "b": b_value - 0.5 * h} + ) + max_index = min(sol_plus.y.shape[1], sol_neg.y.shape[1]) - 1 + dydb_fd = (sol_plus.y[:, :max_index] - sol_neg.y[:, :max_index]) / h + dydb_fd = dydb_fd.transpose().reshape(-1, 1) + + np.testing.assert_array_almost_equal( + dydb_ida[: (2 * max_index), :], dydb_fd + ) + def test_set_atol(self): model = pybamm.lithium_ion.DFN() geometry = model.default_geometry