compdyn / partmc

Particle-resolved stochastic atmospheric aerosol model
http://lagrange.mechse.illinois.edu/partmc/
GNU General Public License v2.0
28 stars 16 forks source link

Look into fixing the shape and size of the Jacobian during CVODE solving #83

Closed mattldawson closed 4 years ago

mattldawson commented 5 years ago

per @cguzman95 Slack convo: in our Jac call on phlex_solver we have this:

// Reset the Jacobian
   SM_NNZ_S(J) = SM_NNZ_S(md->J_init);
   for (int i=0; i<SM_NNZ_S(J); i++) {
   (SM_DATA_S(J))[i] = (realtype)0.0;
   (SM_INDEXVALS_S(J))[i] = (SM_INDEXVALS_S(md->J_init))[i];
   }
   for (int i=0; i<=SM_NP_S(J); i++) {
     (SM_INDEXPTRS_S(J))[i] = (SM_INDEXPTRS_S(md->J_init))[i];
   }

In cvode code, Just before call our Jac cvode calls SUNMatZero(cvdls_mem->A) reseting the data values to zero (on cvode_direct.c) so (SM_DATA_S(J))[i] = (realtype)0.0; is not necessary in fact since cvode is reseting it Plus the interesting is hidden after our jacobian call just after we prepare our jacobian on partmc, cvode look at the jacobian prepared and do this(sunmatrix_sparse.c):

For large cell values at once (tested with 30*30 cells), it seems is going to case 2 or case 3 (most probably 3), because this part is taking like 80% of the time execution. Also each execution takes the same time so is going to only 1 of this case, and some errors Failed model state update: [spec 0] raises (edited)

cguzman95 commented 5 years ago

This issue is fixed on commit f6a083515d55117249bb46ca03283f29be1c16bc, on phlex_solver.c jac_init by modifying this: (SM_INDEXVALS_S(M))[i_elem++] = i_row; by this: (SM_INDEXVALS_S(M))[i_elem++] = i_row+n_dep_var*i_cell;

mattldawson commented 5 years ago

let's still look into avoiding the re-shaping of the Jacobian that CVODE does during solving so we can avoid the Jacobian reset code

cguzman95 commented 5 years ago

It should be related to the CVodeReinit or similar since the Jacobian fails equal times we call solver_run, meaning the indices are right during CVODE solving. However, I have not found a quick solution

cguzman95 commented 4 years ago

I have tested it again and seems solved, CVODE doesn't resize the Jacobian during solving. I believe this was solved time ago when @mattldawson fixed some jacobian calculations, but it seems we forgot to close the issue.

Tested with: CB05 and mock_monarch_1

If this happens again in any test, reopen the issue (to notice that, the calc_jac computation time should be increased notably)