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.

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

Improve numerical stability when integrating the ODE system #68

Closed
4 of 7 tasks
schneiderfelipe opened this issue Nov 26, 2020 · 1 comment · Fixed by #69
Closed
4 of 7 tasks

Improve numerical stability when integrating the ODE system #68

schneiderfelipe opened this issue Nov 26, 2020 · 1 comment · Fixed by #69
Assignees
Labels
bug Something isn't working enhancement New feature or request help wanted Extra attention is needed
Milestone

Comments

@schneiderfelipe
Copy link
Member

schneiderfelipe commented Nov 26, 2020

I encountered some numerical problems when attempting to solve a large system (observe the initial concentration, 1e-5):

$ overreact data/perez-soto2020/RI/BLYP-D4/def2-TZVP/model.k \
    "Benzaldehyde(dcm): 1e-5" "NButylamine(dcm): 1e-5" --plot

gives rise to (observe the explosion in the benzaldehyde concentration and the negative hemiaminal concentration, among other discrepancies):

  no   compound            t = 0 s   t = 1e+11 s
 ────────────────────────────────────────────────
   0   Benzaldehyde(dcm)     0.000     17223.494
   1   NButylamine(dcm)      0.000         0.039
   2   A_N(dcm)              0.000         0.022
   3   A_N_N(dcm)            0.000         0.000
   4   Water(dcm)            0.000        -0.001
   5   A_N_W(dcm)            0.000        -0.000
   6   A_N_N_W(dcm)          0.000        -0.000
   7   A_N_W_W(dcm)          0.000         0.000
   8   TS1_#(dcm)            0.000         0.000
   9   Hemiaminal(dcm)       0.000        -8.834
  10   TS2_#(dcm)            0.000         0.000
  11   I_W(dcm)              0.000        -0.000
  12   TS1N_#(dcm)           0.000         0.000
  13   Int_N(dcm)            0.000         0.000
  14   TS2N_#(dcm)           0.000         0.000
  15   I_N_W(dcm)            0.000        -0.000
  16   TS1W_#(dcm)           0.000         0.000
  17   Int_W(dcm)            0.000        -0.000
  18   TS2W_#(dcm)           0.000         0.000
  19   I_W_W(dcm)            0.000         0.000
  20   TS1WW_#(dcm)          0.000         0.000
  21   Int_W_W(dcm)          0.000         0.000
  22   TS2WW_#(dcm)          0.000         0.000
  23   I_W_W_W(dcm)          0.000        -0.000
  24   Imine(dcm)            0.000         1.273
  25   I_N_W_W(dcm)          0.000         0.000

and the following profile:

explosion

When using 1e-6 as concentration, reaction time becomes sane (0.0001 s versus 1e+11 s), plus the graph becomes flat (no reaction, that is; final concentrations are precisely the same as in the beginning).

When using 1e-4 as concentration, I get an overflow:

/home/schneider/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:336: RuntimeWarning: overflow encountered in multiply
  new_factor = NUM_JAC_FACTOR_INCREASE * factor[ind]
/home/schneider/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:358: RuntimeWarning: overflow encountered in multiply
  factor[max_diff < NUM_JAC_DIFF_SMALL * scale] *= NUM_JAC_FACTOR_INCREASE

Discussion

The above is a problem with the Jacobian, which we currently numerically approximate (EDIT: more precisely, the initial rate is probably poorly estimated). That is not ideal, so I propose a couple of things to mitigate this problem:

  • Allow the user to choose among ODE solvers (see method in scipy.integrate.solve_ivp).
  • Allow the user to select an initial time step.
  • Calculate an analytical Jacobian and use it.
    • Print its eigenvalues at the beginning of the output.

Bear in mind that this is probably related to #65.

EDIT: I think we can consider two other things:

  • the relative tolerance (rtol, a good starting point is 1e-6)
  • the absolute tolerance (atol, a good starting point is 1e-11)

EDIT: Other things to consider:

  • Ensure the diffusion limit is enforced in all cases for reactions in solution (Collins-Kimball theory). This gives an upper bound to the reaction rates, which is the reason this is considered here.
@schneiderfelipe schneiderfelipe added bug Something isn't working enhancement New feature or request help wanted Extra attention is needed labels Nov 26, 2020
@schneiderfelipe schneiderfelipe added this to the 1.0 milestone Nov 26, 2020
@schneiderfelipe schneiderfelipe added this to To do in Chemical kinetics via automation Nov 26, 2020
@schneiderfelipe schneiderfelipe added this to To do in Quality control via automation Nov 26, 2020
@schneiderfelipe schneiderfelipe self-assigned this Nov 26, 2020
@schneiderfelipe schneiderfelipe added this to To do in First publication via automation Nov 26, 2020
@schneiderfelipe schneiderfelipe pinned this issue Nov 26, 2020
@schneiderfelipe
Copy link
Member Author

The most important aspects of this were already included in #69, so I will close this and open specific, lower priority issues.

First publication automation moved this from To do to Done Nov 29, 2020
Quality control automation moved this from To do to Done Nov 29, 2020
Chemical kinetics automation moved this from To do to Done Nov 29, 2020
@schneiderfelipe schneiderfelipe unpinned this issue Nov 29, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working enhancement New feature or request help wanted Extra attention is needed
Projects
Development

Successfully merging a pull request may close this issue.

1 participant