astro-turing / Integrating-diagenetic-equations-using-Python

Reactive-transport model simulating formation of limestone-marl alternations
Apache License 2.0
0 stars 1 forks source link

Merge all branches that use the full power of solve_ivp, i.e. without the py-pde wrapper #4

Closed HannoSpreeuw closed 2 years ago

HannoSpreeuw commented 2 years ago

Currently there are five branches in this repo, including main.

The Fiadeiro-Veronis branch achieves the desired result: it integrates smoothly, without setting any limits on the states. And the end result, after integrating over Tstar, is very similar to Fig 3e from l'Heureux.

It turns out that an explicit (in time) solver (RK23, i.e. Runge Kutta) from _solveivp performs equally well as an implicit (in time) solver like BDF from _solveivp, in the sense that both algorithms produce virtually the same plot after integrating over Tstar; the differences are really minor and can only be seen by blinking the two plots one after the other. This means the implicit solvers from solve_ivp and the Jacobian do not seem to be needed for computing Scenario A from l'Heureux.

Using RK23 it takes 27 minutes on a AMD EPYC 7402P 24-Core Processor, with 48 logical cores to cover Tstar, while with BDF it takes 3 hours and 3 minutes. Towards the end of the integration the memory used is about 12 GB voor RK23 and something similar for BDF. So really huge.

We want to retain the full power of _solveivp and the Jacobian computation for integrating other scenarios that may need implicit (in time) solvers. This means we want to keep a branch that deploys the full power of _solveivp, i.e. _solveivp without the py-pde wrapper, the _Use_solve_ivp_without_py-pdewrapper branch. So we want to merge all branches derived from _Use_solve_ivp_without_py-pdewrapper into _Use_solve_ivp_without_py-pdewrapper.

Yet it would also be good to have a branch that integrates Scenario A fast and with a low memory footprint. This should become the main branch. That fast and efficient integrator for Scenario A can possibly be provided by the py-pde solvers, i.e. using the full py-pde framework, part of which was abandoned in the _Use_solve_ivp_without_py-pdewrapper branch.