Use solve_ivp without the py-pde wrapper, but retain fields and grids from py-pde, for more readable code #44

HannoSpreeuw · 2024-09-03T18:31:53Z

Given the arguments from this closed discussion, I think we should proceed with the methodology from the Use_solve_ivp_without_py-pde_wrapper_branch and merge it into main to include this into our release.

In summary, the Use solve ivp without py pde wrapper branch offers these features:

Faster integrations than main branch.
Allows for all features of solve_ivp, such as a Jacobian sparsity matrix
Event tracking has been implemented, for a number of fields and derived quantities, i.e. the times of sign changes are recorded.

while the unique features of the main branch are less important to us.

…eureux. Expressions become even more involved, unfortunately.

… is the opposite of what is expected. Also end up with equally unstable solutions, i.e. the well-known oscillations along the depth axis. Apparently the sparsity matrix does result in trying to find a better solution, but ultimately fails. Also fixed a bug in Derive_Jacobian.py which luckily was not present in the actual rhs computations. Implemented first steps in setting up the structure for Jacobian computations, so beyond the sparsity matrix. A lot of Numba obviously.

…ent yet. However, it is possible to integrate over longer times without an overflow error. But also running into memory overflows.

…r all Jacobian indices. I had to bypass a couple of Numba's limitations to get that working. However, it ultimately leads to a memory overflow for long integrations. Removing all njit decorators in Compute_jacobian.py makes it faster! So weird. So something is still wrong at the Numba level there. It is also clear now that, even with a Jacobian, the integrations fail and lead to bogus results such as negative cCO3 concentrations and negative values for Phi, so log Phi leads to a "FloatingPointError: invalid value encountered in log". Big bummer. Perhaps try Radau instead of BDF. Or LSODA?

…h slower than the noncompiled version. I still do not understand why and I tried some things to figure this out, e.g. by not using the jac00...jac44 functions, but that did not make a difference. The Jacobian_adjusted_to_reuse_common_factors_and_powers_cleared.txt file should have all the correct Jacobian terms since it comes directly from Derive_Jacobian.py after which I substituted some terms in an editor.

… a vectorize decorator to the overloaded np.heaviside function with the help of numba/numba@2ae154f

…ange has been replaced by range. This is done to investigate if there is an effect on memory use and run time. Less memory use is expected single threaded and perhaps the run times will not increase since integrations seem to scale very poorly with the number of available cores.

…threaded computations.

…ning that off may lead to more accurate results. Also, 'nogil = True' is no longer needed for single core computations.

… in combination with calling functional Jacobians. This was attributed to repeated Numba compilations, possibly from Numba problems with nested functions. That is why the compute_all_Jacobian_elements nested function has been removed (as a separate function), but slowness persisted.

…obians, since that was a long time ago, but as far as I remember none of these runs were successful. This is confirmed by more recent results when we provided a Jacobian sparsity matrix, see commit 2197188: the diagonals should be banded, i.e. have a width of more than 1 element, in order for the Jacobian sparsity matrix to enhance integration. The same should apply for analytical Jacobians. However, the off-diagional elements of an analytical Jacobian matrix will be very hard to compute. This means that at this point we will stop our efforts on deriving analytical Jacobian matrices and only provide Jacobian sparsity matrices. For some background on why Jacobian matrices are banded, please see the literature, e.g. equation 2.58 (page 28) from Finite Diﬀerence Methods Finite Diﬀerence Method for Differential Equationsby Randall J. Leveque (https://edisciplinas.usp.br/pluginfile.php/41896/mod_resource/content/1/LeVeque%20Finite%20Diff.pdf). This example is about the Jacobian for solving the pde describing the motion of a pendulum with a certain mass at the end of a rigid (but massless) bar.

…in'. This means that parameters will be in a separate file and not in ScenarioA.py. Also, we want to include the marlpde folder.

… in the literature.

…als, i.e. the diagonals will have the width of only one depth node. We now know that due to discretisation, there should also be non-zero elements adjacent to the diagonals, but these are very hard to compute.

…be derived after 'eq' has been defined. 'number_of_depths' --> 'Number_of_depths'.

…d in this module. However, it is now calculated in the parameters module, so we no longer need these imports here.

…__' it will complain that it does not know 'no_t_eval'.

…ons, times and metadata using py-pde's FileStorage class seems cumbersome when no tracking is applied. Therefore I reverted to the regular way of saving Numpy arrays in an hdf5 file. All the metadata, i.e. 'stored_parms', which is a single dict, can be stored as well, except for the Jacobian sparsity matrix, since that will give rise to 'TypeError: Object dtype dtype('O') has no native HDF5 equivalent'. This csr_matrix has to be converted to a ndarray first, I reckon.

…eeded for this branch.

…che=False may be noticeably slower for short runs, but sometimes causes a run to halt from start, when a compiled object is missing.

… slightly different data format that solve_ivp returns, i.e. solution.y contains the solutions for the five fields and the [:, -1] indexing gives the last ones across all depths. integrate_equations now gives six return values instead of five. solve_ivp uses 'first_step' instead of 'dt'.

… Numba-based evaluations of the right-hand sides. The Solver dataclass now provides for that. A conditional has been added to check if the 'jac_sparsity' attribute exists. It will not exist for explicit solvers.

…ield, this is more convenient than a single dimension covering all depths for all fields. Also store in this way, as an hdf5. Have 'integrate_equations' only return the final solution, since we only use that for plotting.

…stored results taking excessive disk space.

…of 6. And only the final solutions, which makes comparison with the ground truth somewhat simpler in terms of indexing.

… formatting of floats is often more readable.

…store them in the hdf5 file, one has to iterate over this list and create a separate dataset for each list item, i.e. for each ndarray.

HannoSpreeuw · 2024-09-05T09:49:30Z

I am currently looking at issue #43 , which is about documentation.

Last bullet: saving U at the bottom, this feature will vanish in the merge.

In main it is currently implemented through

    if tracker_parms["track_U_at_bottom"]:
        data_tracker = DataTracker(eq.track_U_at_bottom, \
                               interval = tracker_parms["data_tracker_interval"])

and

    sol, info = eq.solve(state, **solver_parms, tracker=["progress", \
                         storage.tracker(\
                             tracker_parms["progress_tracker_interval"]),\
                         live_plots, data_tracker])

without the py-pde wrapper of solve_ivp we can no longer apply the DataTracker class and tracking of U at the bottom will have to be implemented differently.

@EmiliaJarochowska Pls let me know if tracking (and saving) U at the bottom has to be readded.

HannoSpreeuw · 2024-09-25T12:36:57Z

@EmiliaJarochowska Pls let me know if tracking (and saving) U at the bottom has to be readded.

After discussion: will not be readded anytime soon, instead a similar feature in rhythmite will be deployed.

…en merged into 'main'. Grammar correction. A constant porosity diffusion coefficient is now in all branches. Functional Jacobians turn out not to applicable for this project, because of the discretization. The use of py-pde is now limited to its ScalarField and CartesianGrid.

EmiliaJarochowska · 2024-09-25T19:15:05Z

FYI I managed to crash it:

/opt/homebrew/lib/python3.11/site-packages/scipy/integrate/_ivp/radau.py:401: RuntimeWarning: underflow encountered in nextafter
  min_step = 10 * np.abs(np.nextafter(t, self.direction * np.inf) - t)
 26%|█████████                          | 25866/100000 [03:56<07:03, 174.89it/s]/opt/homebrew/lib/python3.11/site-packages/scipy/integrate/_ivp/radau.py:118: RuntimeWarning: underflow encountered in divide
  dW_norm = norm(dW / scale)

by setting Phi0 to 0.9, Phi00 to 0.8, and b to 2.66667.
However, this also crashed the main so we knew it was prone to it.
So I suggest leaving it like this.

Otherwise works and is indeed very fast.

But the number of solver options and their implications for the numerical methods strengthens the argument that some additional documentation is needed. I will continue this discussion in the respective issue.

EmiliaJarochowska

Pylint highlights some minor issues here and there but the code works so I'll merge now.

EmiliaJarochowska · 2024-09-25T19:43:19Z

I interrupted that run but here are the errors:

  File "/Users/emilia/Documents/Documents - UU101746/GitHub/Integrating-diagenetic-equations-using-Python/marlpde/Evolve_scenario.py", line 211, in <module>
    Plot_results(*integrate_equations(solver_parms, tracker_parms, pde_parms))
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/emilia/Documents/Documents - UU101746/GitHub/Integrating-diagenetic-equations-using-Python/marlpde/Evolve_scenario.py", line 104, in integrate_equations
    sol = solve_ivp(eq.fun if backend=="numpy" else eq.fun_numba, 
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/scipy/integrate/_ivp/ivp.py", line 591, in solve_ivp
    message = solver.step()
              ^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/scipy/integrate/_ivp/base.py", line 181, in step
    success, message = self._step_impl()
                       ^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/scipy/integrate/_ivp/radau.py", line 504, in _step_impl
    f_new = self.fun(t_new, y_new)
            ^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/scipy/integrate/_ivp/base.py", line 136, in fun
    def fun(t, y):

HannoSpreeuw · 2024-09-26T07:44:58Z

FYI I managed to crash it:
/opt/homebrew/lib/python3.11/site-packages/scipy/integrate/_ivp/radau.py:401: RuntimeWarning: underflow encountered in nextafter
  min_step = 10 * np.abs(np.nextafter(t, self.direction * np.inf) - t)
 26%|█████████                          | 25866/100000 [03:56<07:03, 174.89it/s]/opt/homebrew/lib/python3.11/site-packages/scipy/integrate/_ivp/radau.py:118: RuntimeWarning: underflow encountered in divide
  dW_norm = norm(dW / scale)
by setting Phi0 to 0.9, Phi00 to 0.8, and b to 2.66667. However, this also crashed the main so we knew it was prone to it. So I suggest leaving it like this.

Otherwise works and is indeed very fast.

But the number of solver options and their implications for the numerical methods strengthens the argument that some additional documentation is needed. I will continue this discussion in the respective issue.

Yes, I think we discussed those cases as part of issue #36 .

HannoSpreeuw · 2024-09-26T08:04:10Z

Btw, you shared the warning, was there also a crash?

We do not want marlpde to crash on very small time steps, so that is why we have under="warn"
in np.seterr(divide="raise", over="raise", under="warn", invalid="raise").

EmiliaJarochowska · 2024-09-26T14:22:53Z

I am not sure if it crashed, it got stuck at 26% and didn't move, just displayed the warning. Shall I re-try?

HannoSpreeuw · 2024-09-26T15:18:43Z

Yes, please retry, I'd be surprised if it really halted.

EmiliaJarochowska · 2024-09-27T18:03:45Z

It didn't crash, just stuck at 26%. Been at 26% for > 1 h now without any progress. But again, this is a crash case.

HannoSpreeuw · 2024-09-28T08:56:58Z

Okay, this is something we do not want to solve, let's leave it.

EmiliaJarochowska · 2024-09-28T10:07:13Z

Oh, now I discovered that it actually finished. But with the following:

Status = -1 

Success = False

and

Message from solve_ivp = Required step size is less than spacing between numbers.

but yes, it's not to be fixed.

HannoSpreeuw added 30 commits August 29, 2024 16:44

Adjusted the Jacobian calculation to also include equation 6 from l'H…

85c5a92

…eureux. Expressions become even more involved, unfortunately.

First runs using a Jacobian. Unfortunately not a spectacular improvem…

ee0756c

…ent yet. However, it is possible to integrate over longer times without an overflow error. But also running into memory overflows.

Finally managed to jit compile the entire Jacobian function by adding…

783b81c

… a vectorize decorator to the overloaded np.heaviside function with the help of numba/numba@2ae154f

Experimenting with speedups from fastmath = True and pip install icc-rt

0db41b7

'parallel = True' does not make sense in Numba decorators for single …

4e0a266

…threaded computations.

fastmath may lead to an accumulation of errors in integrating, so tur…

fb795b3

…ning that off may lead to more accurate results. Also, 'nogil = True' is no longer needed for single core computations.

Hoped to achieve a speedup with cache=True, but this is not the case.

5e58f4e

We are updating this branch to provide for the same modularity as 'ma…

7f8cde8

…in'. This means that parameters will be in a separate file and not in ScenarioA.py. Also, we want to include the marlpde folder.

A ChatGPT comment is irreproducible, so I looked for a similar remark…

89d285f

… in the literature.

Apparently, one has to use 'default_factory' in this case.

2c14ddf

Compute_Jacobian will yield a Jacobian with only values on the diagon…

0edaf5d

…als, i.e. the diagonals will have the width of only one depth node. We now know that due to discretisation, there should also be non-zero elements adjacent to the diagonals, but these are very hard to compute.

'y0' is derived using a call to a method from 'eq', so 'y0' needs to …

d83275d

…be derived after 'eq' has been defined. 'number_of_depths' --> 'Number_of_depths'.

PyLint seems to prefer f-string formatting here.

23173aa

This variable was too long and therefore renamed.

499c279

This variable needs to be defined before the context can be created.

cd2965e

Again, PyLint seems to prefer f-string formatting here.

24af80f

This branch does not use PDEBase.

d421b08

This branch does not use FieldBase.

e848928

We needed these imports when the Jacobian sparsity matrix was compute…

4172896

…d in this module. However, it is now calculated in the parameters module, so we no longer need these imports here.

PyLint apparently wants this import at the top, i.e. first.

1e73015

PyLint demands imports to be grouped.

8bf3b05

I guess a minimal doctring should be helpful here.

17eeef1

Trailing whitespace.

27eef33

If 't_eval' is defined in the regular way, i.e. not in a '__post_init…

6d8017f

…__' it will complain that it does not know 'no_t_eval'.

HannoSpreeuw added 8 commits September 2, 2024 15:41

ipython is only useful for debugging and prototyping. pandas is not n…

0445519

…eeded for this branch.

Profiling has been done, so the line profiler is no longer needed. ca…

882d9c5

…che=False may be noticeably slower for short runs, but sometimes causes a run to halt from start, when a compiled object is missing.

Clarification for users. Set default 'no_t_eval' really low to avoid …

c339393

…stored results taking excessive disk space.

Now we are back to 'integrate_equations' returning 5 objects instead …

366b523

…of 6. And only the final solutions, which makes comparison with the ground truth somewhat simpler in terms of indexing.

HannoSpreeuw requested a review from EmiliaJarochowska September 3, 2024 18:31

HannoSpreeuw self-assigned this Sep 3, 2024

HannoSpreeuw changed the title ~~Use solve ivp without py pde wrapper~~ Use solve_ivp without the py-pde wrapper, but retain fields and grids from py-pde, for more readable code Sep 4, 2024

HannoSpreeuw and others added 3 commits September 4, 2024 17:50

Store the events as well, it will take little disk space. Exponential…

f5aafc6

… formatting of floats is often more readable.

'sol.t_events' is a list of ndarrays of different sizes. In order to …

477d7e9

…store them in the hdf5 file, one has to iterate over this list and create a separate dataset for each list item, i.e. for each ndarray.

Merge branch 'main' into Use_solve_ivp_without_py-pde_wrapper

b79e81b

HannoSpreeuw mentioned this pull request Sep 25, 2024

Generate documentation #43

Open

3 tasks

EmiliaJarochowska approved these changes Sep 25, 2024

View reviewed changes

EmiliaJarochowska merged commit 0c65692 into main Sep 25, 2024
1 check passed

HannoSpreeuw deleted the Use_solve_ivp_without_py-pde_wrapper branch September 26, 2024 07:45

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Use solve_ivp without the py-pde wrapper, but retain fields and grids from py-pde, for more readable code #44

Use solve_ivp without the py-pde wrapper, but retain fields and grids from py-pde, for more readable code #44

HannoSpreeuw commented Sep 3, 2024

HannoSpreeuw commented Sep 5, 2024

HannoSpreeuw commented Sep 25, 2024

EmiliaJarochowska commented Sep 25, 2024

EmiliaJarochowska left a comment

EmiliaJarochowska commented Sep 25, 2024

HannoSpreeuw commented Sep 26, 2024

HannoSpreeuw commented Sep 26, 2024

EmiliaJarochowska commented Sep 26, 2024

HannoSpreeuw commented Sep 26, 2024

EmiliaJarochowska commented Sep 27, 2024

HannoSpreeuw commented Sep 28, 2024

EmiliaJarochowska commented Sep 28, 2024

Use solve_ivp without the py-pde wrapper, but retain fields and grids from py-pde, for more readable code #44

Use solve_ivp without the py-pde wrapper, but retain fields and grids from py-pde, for more readable code #44

Conversation

HannoSpreeuw commented Sep 3, 2024

HannoSpreeuw commented Sep 5, 2024

HannoSpreeuw commented Sep 25, 2024

EmiliaJarochowska commented Sep 25, 2024

EmiliaJarochowska left a comment

Choose a reason for hiding this comment

EmiliaJarochowska commented Sep 25, 2024

HannoSpreeuw commented Sep 26, 2024

HannoSpreeuw commented Sep 26, 2024

EmiliaJarochowska commented Sep 26, 2024

HannoSpreeuw commented Sep 26, 2024

EmiliaJarochowska commented Sep 27, 2024

HannoSpreeuw commented Sep 28, 2024

EmiliaJarochowska commented Sep 28, 2024