How to continue with the Use_solve_ivp_without_py-pde_wrapper branch?

[HannoSpreeuw]

The Use_solve_ivp_without_py-pde_wrapper branch was created to make most use of the features of solve_ivp.
The py-pde wrapper around solve_ivp limits the use of these features, e.g. one cannot use a Jacobian in a functional form.

The creation of this branch seemed worthwhile when the integration of the diagenetic equations derailed way before T*; i.e. the integrations seemed stiff and the use of implicit solvers with Jacobians in functional form seemed appropriate.

At a later stage in the project a wrong value of b turned out to be detrimental for the stability of the integrations and much more important than the use of implicit solvers: we can now integrate the diagenetic equations for very long timespans without problems using explicit solvers (main branch). The Fiadeiro-Veronis scheme provides for additional stability.

So the question is: do we still need a branch that provides for implicit solvers, i.e. solvers for very stiff equations?

If yes, we need to get this branch at the same level as the main branch, which would involve most of the code enhancements from PR #15.
And possibly a constant porosity diffusion coefficient, which would have an effect on the Jacobian computation as well, since the Jacobian has been derived assuming a varying porosity diffusion coefficient.

[EmiliaJarochowska]

We have good reason to expect that implicit solvers would perform better in our case. But updating this branch is also a substantial work commitment so to make the decision we might first investigate the parts of the system that contribute to stiffness. For example, #24 , #21 and #23 , all related to the parameter b, compressibility and hydraulic conductivity. It might be useful to compare the absolute values with those obtained in Fortran. And they, to some extent, can be compared with empirical values.

There is still some literature work on my side to be done for this. Permeability and its ~linear function hydraulic conductivity can be measured in sediment and thus give us a way to assess whether we get "physical" values. And they are also the main culprits affecting stability.

@HannoSpreeuw mentioned the high memory consumption of solve_ivp. My understanding is that the causes can be investigated, i.e. what part of the integration requires excess steps. I do not know at this point how/if this can be done but perhaps this could help us at a later stage.

[HannoSpreeuw]

I have been wondering why this branch is so slow compared to the main branch, so I profiled it using cProfile with visualisation using SnakeViz and found out that most of the time spent is in Numba routines, which should not be the case; Numba should compile the python code for evaluating the rhs once - i.e. shortly after launching the code - and not every time for each call to the rhs function.
Adding cache=True to the njit decorator alleviates the problem somewhat, but it is still slow.

The repeated Numba compilation does not necessarily explain the memory problem though, but there might be a connection. But no clue about such a connection at this point, I must admit.

I will pursue this investigation to get to the cause of the repeated Numba compilation and fix it. Let's then check if that fix has any effect on memory usage.

[jhidding]

I could help here if you like. Some general remarks:

• Avoid allocating memory inside njit routines
• Are you using the Jacobian? It doesn't seem to be decorated.

[EmiliaJarochowska]

This could be a longer use of the resources of maybe both of you. Perhaps hold on for a few days. I am going to ask all of you for input on a plan of what needs to be done for a manuscript, i.e. which runs remain to be done and what steps remain to release the software with the submission. Perhaps we decide that we don't need to pursue this branch anymore?
I would like to make sure we use your time wisely.

[HannoSpreeuw]

The reason for doing the investigation was marlpde bombing out after some 133978 years. This happens when running from the main branch or from branches close to the main branch. That is an issue we all want to fix, I presume, since we want to know if the oscillations shown by rhythmite are real - and not a numerical artefact. They occur around that time after a relatively fast change in calcite composition, as shown by @csummers25.

The advantage of the Use_solve_ivp_without_py-pde_wrapper branch is that it can handle stiff equations with the use of a Jacobian in a functional form. So if it were fast enough, we may be able to integrate those stiff equations beyond 134,000 years. This would be without clamping since marlpde does not have clamping. So any result beyond 134,000 years is interesting, because we want to know if the rhythmite oscillations are an artefact of clamping.

Another nice feature of this branch is that has all the nice event tracking for sign changes of CA, CC, Phi U and W, so we can get useful information at the time of any fast increase or decrease of CA, i,e. near 134,000 years.

But achieving those insights requires fixing the slowness of this branch.

[HannoSpreeuw]

@jhidding Thanks for those comments. I am going to look into memory allocation, but that is done inside the much faster main branch as well, e.g. rate = np.empty_like(state_data) is a memory allocation inside the jitted pde_rhs in the main branch. (rate is the array that contains the rhs values)

My today's comments relate to runs with method='RK23' in solve_ivp, so that is an explicit solver, i.e. without using a Jacobian.

eq.jac calls Compute_jacobian.Jacobian, which is jitted.

[HannoSpreeuw]

I think I fixed the main issues.

When explicit solvers are used, this branch is about as fast as the main branch. Slightly slower, actually.

Implicit solvers are blazingly fast - even faster then explicit solvers - if and only if a numerically approximated Jacobian is used. You can try "BDF", "Radau" or "LSODA".

Beware that this branch does not have a fixed porosity diffusion coefficient, but we know from past runs that this does not make a huge difference.

[EmiliaJarochowska]

we want to know if the oscillations shown by rhythmite are real - and not a numerical artefact

Isn't it correct that we know by now? I am not sure if we can classify them as "real" or "artifacts". They are produced by the model, but their existence depends on the size of the system and $d_\Phi$, which have nothing to do with conditions in real life. The oscillations might be a true solution to that system, but that system includes a term that has no physical representation and was supposed to only aid the numerics, but turned out the be the key parameter.

Based on that I do not agree with:

That is an issue we all want to fix, I presume, since we want to know if the oscillations shown by rhythmite are real - and not a numerical artefact.

[EmiliaJarochowska]

No. We have a scientific question, which is: can these oscillations be preserved in the geological record and be confounded with climatic cycles. We don't need to know the solutions to this system of equations if the answer to that question is no. It becomes a random set of equations at that point.

[HannoSpreeuw]

I see, "which have nothing to do with conditions in real life" makes it a redundant investigation.

[EmiliaJarochowska]

Sorry I should have made it clear after Wednesday but I was too flummoxed. At this point we just need to prepare the software for publication, I will send an email about it, but it's hard to write it in the middle of an intensive workshop. So it might take up to 2 more days.

[HannoSpreeuw]

For what it's worth, implicit solvers give flat solutions, no osillations.

End of investigation.

[EmiliaJarochowska]

Thanks. We will report it in the paper. And maybe we even need to document/include runs as well. So it is an important piece of information.

[HannoSpreeuw]

Use_solve_ivp_without_py-pde_wrapper has been rebased on main thus offering the same modularity as per commit 366b523.

main and Use_solve_ivp_without_py-pde_wrapper have their own unique features, which cannot be merged, so we will retain them in the upcoming release.

[HannoSpreeuw]

With hindsight, I would like to reconsider....

Given the features of the Use_solve_ivp_without_py-pde_wrapper branch:

• 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.

vs. the features of the main branch:

• Allows for live plotting
• Offers a number of py-pde specific solvers
• Tracking and storage classes provided

then I am in favour of merging Use_solve_ivp_without_py-pde_wrapper into main, since the main branch features listed here are not really important to us.

For me, the most important argument is that we can use all features of solve_ivp, which is not going to happen as long as it is wrapped, given this discussion.

I still appreciate the other classes that py-pde provides, i.e. the fields and grids classes, which allow for clean code. These are used throughout the Use_solve_ivp_without_py-pde_wrapper branch.