From 128e02ea3df36e98f2ebb532a32f4af5ec01dd52 Mon Sep 17 00:00:00 2001 From: martinjrobins Date: Tue, 4 Oct 2022 09:19:35 +0100 Subject: [PATCH 1/6] #2331 fix bug in size of yS numpy array when simulation stops early --- pybamm/solvers/c_solvers/idaklu_casadi.cpp | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/pybamm/solvers/c_solvers/idaklu_casadi.cpp b/pybamm/solvers/c_solvers/idaklu_casadi.cpp index 0f123df4a0..e44ce06d4b 100644 --- a/pybamm/solvers/c_solvers/idaklu_casadi.cpp +++ b/pybamm/solvers/c_solvers/idaklu_casadi.cpp @@ -373,6 +373,14 @@ int sensitivities_casadi(int Ns, realtype t, N_Vector yy, N_Vector yp, // resvalsS now has (∂F/∂p i ) 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; + //} + for (int i = 0; i < np; i++) { //std::cout << "dF/dp = [" << i << "] = ["; //for (int j = 0; j < p_python_functions->number_of_states; j++) { @@ -457,7 +465,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 +483,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 +494,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]); } @@ -609,7 +616,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 +656,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 +672,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 ); From 7f08623457b8a2b5d06f811f7c51fa9d79821454 Mon Sep 17 00:00:00 2001 From: martinjrobins Date: Thu, 6 Oct 2022 15:54:14 +0100 Subject: [PATCH 2/6] #2331 add test for sens with events, fix bug in python-based idaklu solver as well --- pybamm/solvers/c_solvers/idaklu.cpp | 13 +- pybamm/solvers/c_solvers/idaklu_casadi.cpp | 165 +----------------- tests/unit/test_solvers/test_idaklu_solver.py | 76 ++++++++ 3 files changed, 83 insertions(+), 171 deletions(-) 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 e44ce06d4b..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; @@ -373,22 +279,7 @@ int sensitivities_casadi(int Ns, realtype t, N_Vector yy, N_Vector yp, // resvalsS now has (∂F/∂p i ) 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; - //} - 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; @@ -398,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])); @@ -412,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; @@ -554,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; @@ -736,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..e8446fc2db 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -246,6 +246,82 @@ 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 From fe0fb8b62774751aa52464e3a9c5cf078f4a049f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 6 Oct 2022 14:55:29 +0000 Subject: [PATCH 3/6] style: pre-commit fixes --- tests/unit/test_solvers/test_idaklu_solver.py | 34 +++++++++++-------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index e8446fc2db..ba4c810c58 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -282,7 +282,7 @@ def test_sensitivities_with_events(self): model, t_eval, inputs={"a": a_value, "b": b_value}, - calculate_sensitivities=True + calculate_sensitivities=True, ) # test that y[1] remains constant @@ -298,29 +298,33 @@ def test_sensitivities_with_events(self): # 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 - }) + 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) + 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 - }) + 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) + np.testing.assert_array_almost_equal( + dydb_ida[: (2 * max_index), :], dydb_fd + ) def test_set_atol(self): model = pybamm.lithium_ion.DFN() From a6bb1151a3251d8cc972617f8f42e7a49d518986 Mon Sep 17 00:00:00 2001 From: martinjrobins Date: Thu, 6 Oct 2022 17:40:36 +0100 Subject: [PATCH 4/6] #2331 add changlog entry --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8c2c65ce0b..0d5e5330ff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,7 +9,7 @@ - For experiments, the simulation now automatically checks and skips steps that cannot be performed (e.g. "Charge at 1C until 4.2V" from 100% SOC) ([#2212](https://github.com/pybamm-team/PyBaMM/pull/2212)) ## 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)) - Arrhenius function for `nmc_OKane2022` positive electrode actually gets used now ([#2309](https://github.com/pybamm-team/PyBaMM/pull/2309)) - Added `SEI on cracks` to loop over all interfacial reactions ([#2262](https://github.com/pybamm-team/PyBaMM/pull/2262)) - Fixed `X-averaged SEI on cracks concentration` so it's an average over x only, not y and z ([#2262](https://github.com/pybamm-team/PyBaMM/pull/2262)) From 25d1d0ba67b442bffcb88e61f1acb58a9691fd08 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 6 Oct 2022 16:43:38 +0000 Subject: [PATCH 5/6] style: pre-commit fixes --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0d5e5330ff..cf20c50f1e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ - For experiments, the simulation now automatically checks and skips steps that cannot be performed (e.g. "Charge at 1C until 4.2V" from 100% SOC) ([#2212](https://github.com/pybamm-team/PyBaMM/pull/2212)) ## 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)) - Arrhenius function for `nmc_OKane2022` positive electrode actually gets used now ([#2309](https://github.com/pybamm-team/PyBaMM/pull/2309)) - Added `SEI on cracks` to loop over all interfacial reactions ([#2262](https://github.com/pybamm-team/PyBaMM/pull/2262)) From 749019d1b012109336ab42aba9dd96bd5687a3fb Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 6 Oct 2022 17:32:55 -0400 Subject: [PATCH 6/6] Update CHANGELOG.md --- CHANGELOG.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cf20c50f1e..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 @@ -10,7 +14,6 @@ ## 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)) - Arrhenius function for `nmc_OKane2022` positive electrode actually gets used now ([#2309](https://github.com/pybamm-team/PyBaMM/pull/2309)) - Added `SEI on cracks` to loop over all interfacial reactions ([#2262](https://github.com/pybamm-team/PyBaMM/pull/2262)) - Fixed `X-averaged SEI on cracks concentration` so it's an average over x only, not y and z ([#2262](https://github.com/pybamm-team/PyBaMM/pull/2262))