Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

#2331 fix bug in size of yS numpy array when simulation stops early #2337

Merged
merged 6 commits into from
Oct 6, 2022
Merged

#2331 fix bug in size of yS numpy array when simulation stops early #2337

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
128e02e
#2331 fix bug in size of yS numpy array when simulation stops early
martinjrobins Oct 4, 2022
7f08623
#2331 add test for sens with events, fix bug in python-based idaklu s…
martinjrobins Oct 6, 2022
fe0fb8b
style: pre-commit fixes
pre-commit-ci[bot] Oct 6, 2022
a6bb115
#2331 add changlog entry
martinjrobins Oct 6, 2022
25d1d0b
style: pre-commit fixes
pre-commit-ci[bot] Oct 6, 2022
749019d
Update CHANGELOG.md
valentinsulzer Oct 6, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 6 additions & 7 deletions pybamm/solvers/c_solvers/idaklu.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<realtype *> ySval(number_of_parameters);
int retval;
SUNMatrix J;
SUNLinearSolver LS;
Expand All @@ -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;
Expand All @@ -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]);
}
Expand Down Expand Up @@ -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];
}
}

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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<ptrdiff_t>{number_of_parameters, t_i, number_of_states},
std::vector<ptrdiff_t>{number_of_parameters, number_of_timesteps, number_of_states},
&yS_return[0]
);

Expand Down
170 changes: 7 additions & 163 deletions pybamm/solvers/c_solvers/idaklu_casadi.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -104,59 +95,25 @@ int residual_casadi(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr,
CasadiFunctions *p_python_functions =
static_cast<CasadiFunctions *>(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);
p_python_functions->rhs_alg.m_arg[2] = static_cast<realtype *>(input_buf.ptr);
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;
}

Expand Down Expand Up @@ -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];
}

Expand All @@ -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];
}

Expand All @@ -278,17 +227,6 @@ int events_casadi(realtype t, N_Vector yy, N_Vector yp, realtype *events_ptr,
CasadiFunctions *p_python_functions =
static_cast<CasadiFunctions *>(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;
Expand All @@ -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);
}

Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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]));

Expand All @@ -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;
Expand Down Expand Up @@ -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<realtype *> ySval(number_of_parameters);
int retval;
SUNMatrix J;
SUNLinearSolver LS;
Expand All @@ -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;
Expand All @@ -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]);
}
Expand Down Expand Up @@ -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<realtype> 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<double> 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;
Expand Down Expand Up @@ -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];
}
}

Expand Down Expand Up @@ -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;
Expand All @@ -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<ptrdiff_t>{number_of_parameters, t_i, number_of_states},
std::vector<ptrdiff_t>{number_of_parameters, number_of_timesteps, number_of_states},
&yS_return[0], free_yS_when_done
);

Expand Down Expand Up @@ -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;
}

Loading