Bit of 'digging' into the stratigraphy

[elbeejay]

This is not meant to ever be merged, but rather to prompt a discussion about what we should do regarding stratigraphy recording (this grew out of efforts to close #43).

'old' method

The original DeltaRCM stratigraphy process involved initializing a 3D array and populating it over time with sand fraction values using a 2D array sand_frac and the vertical locations to be recorded were assigned using another 2D array topz. The stratigraphy array was updated every timestep after sediment routing was completed. This implementation was carried over to relauzon's python DeltaRCM which was the genesis for what is now pyDeltaRCM_WMT.

'new' method

I suspect, that when the code was being updated and improved, this initialization of a full 3D array and writing to it every timestep was identified as a real processing bottleneck. Which is why the current method carefully initializes smaller sparse arrays to hold the strata sand content and depth that are only written to and extended at timesteps that are output (intervals of save_dt). In this 'new' method, the sand fractions and depth of sand is maintained throughout the simulation, and a 'final' stratigraphy can be back calculated using strat_preprocess.py (added in this PR as tools.strata_postprocess) to cross-cut between the sand fraction arrays and the strata depth arrays to get at a final saved stratigraphy.

comparing the two

This jump from the original implementation to the "new" or "sparse" one does save a lot of computational time. But, how do we know if it is accurate? That is why I put together this PR to resurrect the 'old' method and directly compare it to the new one. A direct comparison of the stratigraphy output via these two methods shows us that that are not in-fact identical (image). Fortunately, the strata bits of code are very much independent from the rest of the model, so this inconsistency can be explored in isolation (topography and discharge remain the same).

So this is kinda the status of how the model handles stratigraphy. I am still thinking through what the 'best' path forward needs to be so that we are as confident as we can be in the pyDeltaRCM stratigraphy that is modeled. Suggestions, comments and ideas are appreciated.