diff --git a/CHANGELOG.md b/CHANGELOG.md index 18bf754946..cef03a025e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ ## 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)) ## Bug Fixes diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp index 36c2872c3e..92eede3643 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, yp, 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 *ypS; // cppcheck-suppress unusedStructMember N_Vector id; // rhs_alg_id @@ -70,6 +70,7 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver vector res_dvar_dp; bool const sensitivity; // cppcheck-suppress unusedStructMember bool const save_outputs_only; // 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 @@ -158,6 +159,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 313c4ce12a..56f546facf 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl @@ -83,6 +83,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 @@ -92,12 +96,14 @@ void IDAKLUSolverOpenMP::AllocateVectors() { if (setup_opts.num_threads == 1) { yy = N_VNew_Serial(number_of_states, sunctx); yp = 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); yp = 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); } @@ -289,9 +295,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) @@ -312,6 +322,7 @@ IDAKLUSolverOpenMP::~IDAKLUSolverOpenMP() { N_VDestroy(avtol); N_VDestroy(yy); N_VDestroy(yp); + N_VDestroy(y_cache); N_VDestroy(id); if (sensitivity) { @@ -386,21 +397,16 @@ SolutionData IDAKLUSolverOpenMP::solve( SetSolverOptions(); - CheckErrors(IDAReInit(ida_mem, t0, yy, yp)); - if (sensitivity) { - CheckErrors(IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, ypS)); - } - // 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); } if (sensitivity) { @@ -480,12 +486,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, yp)); - if (sensitivity) { - CheckErrors(IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, ypS)); - } + ReinitializeIntegrator(t_val); + ConsistentInitialization(t_val, t_eval_next, IDA_YA_YDP_INIT); } t_prev = t_val; @@ -563,6 +565,53 @@ void IDAKLUSolverOpenMP::ExtendAdaptiveArrays() { } } +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 ea3903b139..32048d89c0 100644 --- a/src/pybamm/solvers/idaklu_solver.py +++ b/src/pybamm/solvers/idaklu_solver.py @@ -988,7 +988,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