Investigate the differences between marlpde and rhytmite when the same solver and initial settings are used

[EmiliaJarochowska]

We want to understand why marlpde and rhytmite give different results for runs longer than 1 dt when the same solver settings are used.

This can be done by fixing the solver and the parameters and focusing on why the results still differ. I suggest focusing on Euler using a fixed time time step of 1e-6. But if needed, a comparison of RK23, RK45 or adaptive could also be done. Important is to do one comparison at a time.

This should be done with the corrected typo as in #39 so using the Oscillations_like_Fortran branch. Additional adjustment can be made, either to marlpde or to rhytmite, e.g. the clipping of values, to pinpoint the causes of differences.

Here is a comparison using the Euler solver after 100 dt, made by Charlotte (email from 13 June):

[EmiliaJarochowska]

The differences are small but you can use the code for plotting residues in https://github.com/MindTheGap-ERC/Cross-comparison to investigate them

[HannoSpreeuw]

Okay, I guess this committed comparison plot, from the day before, i.e. June 12 was also made using the main marlpde branch, i.e. with the correct formula for the porosity diffusion coefficient instead of the restored typo.

Including the default solver settings from the main marlpde branch.

I will use new folders to upload comparisons coming from anything different from the main marlpde branch.

[HannoSpreeuw]

Better resemblance now near the surface and away from the ADZ.

See these two commits for some background.

[HannoSpreeuw]

This is also informative.

[EmiliaJarochowska]

That is indeed very good, considering that we are interested in what is exported from the system, some discrepancy in the ADZ is probably acceptable, as long as the values converge in the lower part of the system. Any ideas what the causes of the remaining differences might be?

[HannoSpreeuw]

I guess the remaining differences can be explained by either:

• slightly different position of the rhythmite ADZ wrt marlpde, since the rhythmite grid is shifted wrt marlpde: marlpde effectuates boundary conditions at the upper edge of the first cell - at index = -0.5 for the surface boundary condition - and at the lower edge of the last cell while rhythmite does so at the centers of the first and last cells.
• different calculation of the CA and CC spatial derivatives

I can continue investigating this in a week from now.

[HannoSpreeuw]

In order to resolve the remaining differences the simplest approach would be to align the rhythmite and marlpde grids, since they are half a cell apart. Aligning the grids would likely result in a closer resemblance than the linear interpolation added to the comparison code through this commit.

Aligning the grids is impossible without changing the way rhythmite effectuates the boundary conditions, since these are imposed at cell centers, not at cell edges.

[HannoSpreeuw]

The easiest way to "resolve" the differences is to revert to a much higher spatial resolution, e.g. a grid of 1000 instead of 200 depths, and check if the differences in the above figure diminish.

[HannoSpreeuw]

However, nnx >= 400 will make rhythmite exit with nan encountered, exiting.

[HannoSpreeuw]

Tried a different approach. No ADZ ---> No clamping is needed in rhythmite. This makes an integration over T* possible. These are the differences between rhythmite and marlpde.

See this commit for the code.

[EmiliaJarochowska]

Here's what I received from Charlotte:

Rhythmite can be run at higher resolutions as it is, but the timestep needs to be decreased to keep it stable, I found that with dt=10-7, it can do at least 800 nnx without issue.

[HannoSpreeuw]

Okay, thanks, I hadn't changed the time step.

For such a large depth grid one will probably need to run this on a powerful cluster node, if one wants to integrate over an interval appreciably longer than 1e-4 times T*.

[EmiliaJarochowska]

Turns out Charlotte has ran the comparison for the steady-state solution (Fig. 3) with different grid sizes (200 to 1000 by 100) and the agreement between rhytmite and marlpde is excellent. There are only very tiny differences at the surface of the system, which is not a problem.

[EmiliaJarochowska]

@HannoSpreeuw observed that if clamping of non-physical values is removed from rhytmite, then integration over longer times (as in T*) is not possible, unless one also switches off ADZ. So clamping of unphysical values seems crucial, but is it possible to obtain solutions without it? According to Cedric, if intermittent values are unphysical but the final result is, then they are acceptable as an approximation of the solver. I added @csummers25 here btw :-)

[HannoSpreeuw]

Here is yet another way of assessing differences between marlpde and rhythmite. This time a 250 ka run, with initial and surface porosities of 0.8. In marlpde/parameters.py solver: str = "explicit" scheme: str = "rk"

After some 133978 years one runs into

RuntimeError: Time step below 1e-10 54%|███████████████████████████████████████████▍ | 10.15751/18.953752843062926

Now this time of 133978 years may coincide with the calcite bump prior to oscillations from the slides from @csummers25

What does this mean? marlpde tries to enforce a tolerance=0.0001 but that would yield a time step < 1e-10 when the calcite composition becomes almost 1?

But there is a need to verify that also marlpde's calcite composition at the time of bombing out (~133978 years) is almost 1, i.e. in agreement with rhythmite's value at that same time.

I will try an Eulerian scheme as well.

[HannoSpreeuw]

Same result for scheme: str = "euler", slighty earlier, i.e.

54% --> 53%

10.15751 --> 10.12714

[HannoSpreeuw]

With solver: str = "scipy" in marlpde/parameters.py and method="LSODA", first_step=1e-6, atol=1e-3:

RuntimeError: Encountered Not-A-Number (NaN) in evolution 53%|███████████████████████████████████████████▍ | 10.04276/18.953752843062926

So a similar result.

[HannoSpreeuw]

A possible way out would be to add clamping to physically feasible values to marlpde and see if it can integrate beyond this point.

[HannoSpreeuw]

Clamping does not seem to be needed to integrate over 250,000 years if implicit solvers -with numerically approximated instead of functional Jacobians - in marlpde's updated Use_solve_ivp_without_py-pde_wrapper branch is used.

I find these solutions at 50 cm depth over that time interval - i.e. no oscillations:

[HannoSpreeuw]

With a constant porosity diffusion coefficient and the implicit "BDF" solver, those profiles at 50 cm depth are even flatter. Using that solver from solve_ivp this result can be obtained from the main branch as well as from the Use_solve_ivp_without_py-pde_wrapper branch; the latter has been updated to run with a constant porosity diffusion coefficient.

[HannoSpreeuw]

Does rhythmite have a varying porosity diffusion coefficient?

[EmiliaJarochowska]

Yes, please see https://github.com/cedrict/rhythmite/blob/257136bec97acc283dbc9e682c7fcf6b8669b387/DiagenesisModel.py#L142-L144

[HannoSpreeuw]

Almost the same link.

[csummers25]

It doesn't vary during the simulation, the value is just calculated from the other parameters initially

[HannoSpreeuw]

Ah, now I see, sorry for the oversight.

[HannoSpreeuw]

@csummers25 @EmiliaJarochowska

I added some commits to the Oscillations_like_Fortran branch.

Pivotal in making marlpde produce oscillations are two things:

1) Introduce the typo. In marlpde the original typo was mimicked, but the effect should be the same. 2) Make sure that backward spatial differencing for Aragonite and Calcite is replaced by forward differencing to retain an upwind solution when U changes sign. This required reverting commit 0d8191f , such that the "I guess this seems reasonable" boundary condition - second spatial derivative, i.e. curvature = 0 - for aragonite and calcite is reintroduced.

This is the result from commit c1d98ae (second latest commit from the Oscillations_like_Fortran branch.

[HannoSpreeuw]

In an effort to mimick the bottom boundary condition for aragonite and calcite used by rhyhtmite, I added commit 73b91b4.

See the log message of that commit to understand why it may come closer to the rhythmite solution for the missing bottom boundary conditions for aragonite and calcite than the curvature=0 bottom boundary condition from commit c1d98ae.

The result is similar, though.

[HannoSpreeuw]

Plotting script.

import numpy as np
import h5py
import pylab as pyl
import os

Tstar = 13_190
pyl.clf()
hf = h5py.File('LMAHeureuxPorosityDiff.hdf5', 'r')
data = np.array(hf["data"])
times = np.array(hf["times"])
times *= Tstar
ms = 3
pyl.title("Compositions, concentrations and porosity over 250 ka at 50 cm depth, commit hash", fontsize=9)
pyl.plot(times, data[:, 0, 20], "v", ms = ms, label = "CA")
pyl.plot(times, data[:, 1, 20], "^", ms = ms, label = "CC")
pyl.plot(times, data[:, 2, 20],  ">", ms = ms, label = "cCa")
pyl.plot(times, data[:, 3, 20], "<", ms = ms, label = "cCO3" )
pyl.plot(times, data[:, 4, 20], "o", ms = ms, label = "Phi")
pyl.legend(loc='lower left')
pyl.savefig("Oscillations_from_marlpde_using_Oscillations_like_Fortran_branch_commit_hash.png", bbox_inches="tight")

[HannoSpreeuw]

Commit 24cd1b5 has a slope=0 imposed on the aragonite and calcite profiles at the bottom, instead of a second order derivative=0 or instead of slope=adjacent_value, which should mimick rhythmite.

The oscillations are again similar.

[HannoSpreeuw]

Doubling the grid dimension, from 5m to 10m, while retaining the depth resolution makes the oscillations vanish.

[HannoSpreeuw]

And a similar result with a curvature=0 ad-hoc bottom boundary condition for aragonite and calcite:

[EmiliaJarochowska]

Thank you for this thorough investigation. How to apply the conclusions? I would like to suggest reverting the typo into the code of marlpde to the correct form (**3) and adjusting b accordingly, as it has been done in rhytmite, for consistency. As the purpose of marlpde was to reproduce the oscillations, I would like to merge a version of the boundary conditions that allow that. @HannoSpreeuw you have found several, in the pursuit of mimicking rhytmite. Would this commit be the right version of the code to merge?

[HannoSpreeuw]

Seems good to merge with main after reverting the typo etc.

The most suitable commit to merge with main depends on the result of this discussion, I guess.

[HannoSpreeuw]

I am currently running marlpde with a debugger - and breakpoints in pde/grids/boundaries/local.py - to figure out how

        self.bc_CA = [{"value": CA0}, 
                      {"derivative_expression" : \
                       lambda adjacent_value, dx, x, t: adjacent_value}]
        self.bc_CC = [{"value": CC0}, 
                      {"derivative_expression" : \
                       lambda adjacent_value, dx, x, t: adjacent_value}]   

propagates, i.e. how marlpde computes spatial derivatives at the bottom for aragonite and calcite.

The debugger says adjacent_value=array(0.).

Also value_func(adjacent_value, dx, *args) from return dx * value_func(adjacent_value, dx, *args) + adjacent_value is array(0.).

[HannoSpreeuw]

Similar results from

        self.bc_CA = [{"value": CA0}, 
                      {"value_expression" : \
                       lambda adjacent_value, dx, x, t: adjacent_value}]
        self.bc_CC = [{"value": CC0}, 
                      {"value_expression" : \
                       lambda adjacent_value, dx, x, t: adjacent_value}] 

i.e. adjacent_value=array(0.) and value_func(adjacent_value, *args)=array(0.). (value_func(adjacent_value, *args) from this line.)

With "derivative_expression", I could understand why the spatial derivatives at the bottom should be zero, since we started off with flat initial distribiutions for calcite and aragonite, but for "value_expression" I don't understand the zeros, e.g. why adjacent_value should be zero. I would expect adjacent_value=CAIni and adjacent_value=CCIni.

[HannoSpreeuw]

Anyway, I just realised that if we are to mimic rhythmite, which has the spatial derivatives at the deepest node equal to spatial derivatives at the second deepest node, this is effectively a second spatial derivative=0, i.e. we are back at the curvature=0 bottom boundary condition for aragonite and calcite. This is commit c1d98ae.

So I am going for that commit - with some updates - to be merged with main.

[EmiliaJarochowska]

Thanks for checking this. So we can go back to that commit, but please use the correct (as per the article) version of the porosity diffusion coefficient (**3) and adjust b instead

[HannoSpreeuw]

But set a default value of b=5e-4 when merging with main, I reckon?

[EmiliaJarochowska]

As mentioned in the PR, it has no meaning if you don't use the typo, so hardly a "default" but yes. My concern was that if a user wants to reproduce the oscillations, they will need to change it. So perhaps there is some value in helping the users through the maze of the parameters by setting the defaults as used in our runs for oscillations. But I won't die on this hill.

[HannoSpreeuw]

@EmiliaJarochowska I guess this investigation has been completed and we can close this issue.