diff --git a/CHANGELOG.md b/CHANGELOG.md index ebe4abdd66..cef03a025e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,15 +4,17 @@ - Added sensitivity calculation support for `pybamm.Simulation` and `pybamm.Experiment` ([#4415](https://github.com/pybamm-team/PyBaMM/pull/4415)) - Added OpenMP parallelization to IDAKLU solver for lists of input parameters ([#4449](https://github.com/pybamm-team/PyBaMM/pull/4449)) +- Added phase-dependent particle options to LAM #4369 - Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode (`SplitOCVR`). ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330)) ## Optimizations +- Improved performance of initialization and reinitialization of ODEs in the (`IDAKLUSolver`). ([#4453](https://github.com/pybamm-team/PyBaMM/pull/4453)) - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416)) -## Features +## Bug Fixes -- Added phase-dependent particle options to LAM #4369 +- Fixed bug where IDAKLU solver failed when `output variables` were specified and an extrapolation event is present. ([#4440](https://github.com/pybamm-team/PyBaMM/pull/4440)) ## Breaking changes diff --git a/src/pybamm/solvers/base_solver.py b/src/pybamm/solvers/base_solver.py index dc1ac0cd72..6334169cf0 100644 --- a/src/pybamm/solvers/base_solver.py +++ b/src/pybamm/solvers/base_solver.py @@ -1489,8 +1489,12 @@ def check_extrapolation(self, solution, events): # second pass: check if the extrapolation events are within the tolerance last_state = solution.last_state - t = last_state.all_ts[0][0] - y = last_state.all_ys[0][:, 0] + if solution.t_event: + t = solution.t_event[0] + y = solution.y_event[:, 0] + else: + t = last_state.all_ts[0][0] + y = last_state.all_ys[0][:, 0] inputs = last_state.all_inputs[0] if isinstance(y, casadi.DM): diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp index 092d591d9e..7cef900f5c 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp @@ -54,7 +54,7 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver int const number_of_events; // cppcheck-suppress unusedStructMember int number_of_timesteps; int precon_type; // cppcheck-suppress unusedStructMember - N_Vector yy, yyp, avtol; // y, y', and absolute tolerance + N_Vector yy, yp, y_cache, avtol; // y, y', y cache vector, and absolute tolerance N_Vector *yyS; // cppcheck-suppress unusedStructMember N_Vector *yypS; // cppcheck-suppress unusedStructMember N_Vector id; // rhs_alg_id @@ -71,6 +71,7 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver bool const sensitivity; // cppcheck-suppress unusedStructMember bool const save_outputs_only; // cppcheck-suppress unusedStructMember bool save_hermite; // cppcheck-suppress unusedStructMember + bool is_ODE; // cppcheck-suppress unusedStructMember int length_of_return_vector; // cppcheck-suppress unusedStructMember vector t; // cppcheck-suppress unusedStructMember vector> y; // cppcheck-suppress unusedStructMember @@ -166,6 +167,32 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver */ void PrintStats(); + /** + * @brief Set a consistent initialization for ODEs + */ + void ReinitializeIntegrator(const realtype& t_val); + + /** + * @brief Set a consistent initialization for the system of equations + */ + void ConsistentInitialization( + const realtype& t_val, + const realtype& t_next, + const int& icopt); + + /** + * @brief Set a consistent initialization for DAEs + */ + void ConsistentInitializationDAE( + const realtype& t_val, + const realtype& t_next, + const int& icopt); + + /** + * @brief Set a consistent initialization for ODEs + */ + void ConsistentInitializationODE(const realtype& t_val); + /** * @brief Extend the adaptive arrays by 1 */ diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl index 16f64e6df9..ab39144390 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl @@ -86,6 +86,10 @@ IDAKLUSolverOpenMP::IDAKLUSolverOpenMP( if (this->setup_opts.preconditioner != "none") { precon_type = SUN_PREC_LEFT; } + + // The default is to solve a DAE for generality. This may be changed + // to an ODE during the Initialize() call + is_ODE = false; } template @@ -95,12 +99,14 @@ void IDAKLUSolverOpenMP::AllocateVectors() { if (setup_opts.num_threads == 1) { yy = N_VNew_Serial(number_of_states, sunctx); yyp = N_VNew_Serial(number_of_states, sunctx); + y_cache = N_VNew_Serial(number_of_states, sunctx); avtol = N_VNew_Serial(number_of_states, sunctx); id = N_VNew_Serial(number_of_states, sunctx); } else { DEBUG("IDAKLUSolverOpenMP::AllocateVectors OpenMP"); yy = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); yyp = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); + y_cache = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); avtol = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); id = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); } @@ -312,9 +318,13 @@ void IDAKLUSolverOpenMP::Initialize() { realtype *id_val; id_val = N_VGetArrayPointer(id); - int ii; - for (ii = 0; ii < number_of_states; ii++) { + // Determine if the system is an ODE + is_ODE = number_of_states > 0; + for (int ii = 0; ii < number_of_states; ii++) { id_val[ii] = id_np_val[ii]; + // check if id_val[ii] approximately equals 1 (>0.999) handles + // cases where id_val[ii] is not exactly 1 due to numerical errors + is_ODE &= id_val[ii] > 0.999; } // Variable types: differential (1) and algebraic (0) @@ -335,6 +345,7 @@ IDAKLUSolverOpenMP::~IDAKLUSolverOpenMP() { N_VDestroy(avtol); N_VDestroy(yy); N_VDestroy(yyp); + N_VDestroy(y_cache); N_VDestroy(id); if (sensitivity) { @@ -367,9 +378,7 @@ SolutionData IDAKLUSolverOpenMP::solve( !save_outputs_only ); - // if (t.size() < number_of_evals + number_of_interps) { InitializeStorage(number_of_evals + number_of_interps); - // } int i_save = 0; @@ -415,21 +424,16 @@ SolutionData IDAKLUSolverOpenMP::solve( SetSolverOptions(); - CheckErrors(IDAReInit(ida_mem, t0, yy, yyp)); - if (sensitivity) { - CheckErrors(IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, yypS)); - } - // Prepare first time step i_eval = 1; realtype t_eval_next = t_eval[i_eval]; + // Consistent initialization + ReinitializeIntegrator(t0); int const init_type = solver_opts.init_all_y_ic ? IDA_Y_INIT : IDA_YA_YDP_INIT; if (solver_opts.calc_ic) { - DEBUG("IDACalcIC"); - // IDACalcIC will throw a warning if it fails to find initial conditions - IDACalcIC(ida_mem, init_type, t_eval_next); + ConsistentInitialization(t0, t_eval_next, init_type); } // Set the initial stop time @@ -521,12 +525,8 @@ SolutionData IDAKLUSolverOpenMP::solve( CheckErrors(IDASetStopTime(ida_mem, t_eval_next)); // Reinitialize the solver to deal with the discontinuity at t = t_val. - // We must reinitialize the algebraic terms, so do not use init_type. - IDACalcIC(ida_mem, IDA_YA_YDP_INIT, t_eval_next); - CheckErrors(IDAReInit(ida_mem, t_val, yy, yyp)); - if (sensitivity) { - CheckErrors(IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, yypS)); - } + ReinitializeIntegrator(t_val); + ConsistentInitialization(t_val, t_eval_next, IDA_YA_YDP_INIT); } t_prev = t_val; @@ -666,6 +666,53 @@ void IDAKLUSolverOpenMP::ExtendHermiteArrays() { } } +template +void IDAKLUSolverOpenMP::ReinitializeIntegrator(const realtype& t_val) { + DEBUG("IDAKLUSolver::ReinitializeIntegrator"); + CheckErrors(IDAReInit(ida_mem, t_val, yy, yp)); + if (sensitivity) { + CheckErrors(IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, ypS)); + } +} + +template +void IDAKLUSolverOpenMP::ConsistentInitialization( + const realtype& t_val, + const realtype& t_next, + const int& icopt) { + DEBUG("IDAKLUSolver::ConsistentInitialization"); + + if (is_ODE && icopt == IDA_YA_YDP_INIT) { + ConsistentInitializationODE(t_val); + } else { + ConsistentInitializationDAE(t_val, t_next, icopt); + } +} + +template +void IDAKLUSolverOpenMP::ConsistentInitializationDAE( + const realtype& t_val, + const realtype& t_next, + const int& icopt) { + DEBUG("IDAKLUSolver::ConsistentInitializationDAE"); + IDACalcIC(ida_mem, icopt, t_next); +} + +template +void IDAKLUSolverOpenMP::ConsistentInitializationODE( + const realtype& t_val) { + DEBUG("IDAKLUSolver::ConsistentInitializationODE"); + + // For ODEs where the mass matrix M = I, we can simplify the problem + // by analytically computing the yp values. If we take our implicit + // DAE system res(t,y,yp) = f(t,y) - I*yp, then yp = res(t,y,0). This + // avoids an expensive call to IDACalcIC. + realtype *y_cache_val = N_VGetArrayPointer(y_cache); + std::memset(y_cache_val, 0, number_of_states * sizeof(realtype)); + // Overwrite yp + residual_eval(t_val, yy, y_cache, yp, functions.get()); +} + template void IDAKLUSolverOpenMP::SetStep( realtype &tval, diff --git a/src/pybamm/solvers/idaklu_solver.py b/src/pybamm/solvers/idaklu_solver.py index 10d633c3b0..058bfc0a1a 100644 --- a/src/pybamm/solvers/idaklu_solver.py +++ b/src/pybamm/solvers/idaklu_solver.py @@ -1011,7 +1011,7 @@ def _rhs_dot_consistent_initialization(self, y0, model, time, inputs_dict): rhs0 = rhs_alg0[: model.len_rhs] - # for the differential terms, ydot = -M^-1 * (rhs) + # for the differential terms, ydot = M^-1 * (rhs) ydot0[: model.len_rhs] = model.mass_matrix_inv.entries @ rhs0 return ydot0 diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index 0697f64786..f4028215c1 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -1159,3 +1159,24 @@ def test_python_idaklu_deprecation_errors(self): match="Unsupported evaluation engine for convert_to_format=jax", ): _ = solver.solve(model, t_eval) + + def test_extrapolation_events_with_output_variables(self): + # Make sure the extrapolation checks work with output variables + model = pybamm.BaseModel() + v = pybamm.Variable("v") + c = pybamm.Variable("c") + model.variables = {"v": v, "c": c} + model.rhs = {v: -1, c: 0} + model.initial_conditions = {v: 1, c: 2} + model.events.append( + pybamm.Event( + "Triggered event", + v - 0.5, + pybamm.EventType.INTERPOLANT_EXTRAPOLATION, + ) + ) + solver = pybamm.IDAKLUSolver(output_variables=["c"]) + solver.set_up(model) + + with pytest.warns(pybamm.SolverWarning, match="extrapolation occurred for"): + solver.solve(model, t_eval=[0, 1])