• update release notes in the docs!

Changes:

• refactor stream power law implementation (decouple from flow routing) #23
• rename Union internal subroutine (could cause name conflicts when building the library with f2py) #22

New features:

• new routines for computing curvature and slope #15
• new component for marine sediment transport and deposition #20

Bug fixes and improvements:

• fixed VTK files export on Windows #7
• improved efficiency of stream power law erosion and flexure #8
• fixed bug in lake filling algorithm #12
• fixed bug in ADI implementation of diffusion #17
• fixed bug in Strati subroutine #18

There might be an issue introduced between versions 2.7 and 2.8.

Taking the example in the notebook FastScape_test included in the examples, plotting the topography after half the total simulation time,

The results for version 2.7 and 2.8 in MacOS, and 2.7 in Linux:

And the results for version 2.8 in Linux:

The model doesn't crash, but it seems that SPL sediment transport is not working properly, only for Linux and version 2.8.

This is weird! Maybe it has something to do with how the Python bindings are compiled with F2PY. I haven't checked the output from Fortran examples in Linux, though. Maybe related to fastscape-lem/fastscape#7.

Any idea @jeanbraun ?

Hi,
I've been having an issue with the newest build where the VTK files do not show up in paraview on a GNU system. This seems to be an issue with the removal of "convert=big_endian" from the open function that writes the files, as setting this again fixes the issue.

From my understanding, it looks like cmakelists now sets the flag "-fconvert=big-endian", which should have been applied to the function and done the conversion without specifying it, but for some reason it seems this isn't occurring. The flags themselves seem to be included in the build, as currently with "-std=f2008" included adding "convert=big_endian" to the open function results in an error, but if it is removed I can again specify it like in previous builds. Because of this I'm not sure exactly why the convert flag isn't working, any ideas what the issue may be?

Thank you for your time!

According to the online documentation for the slope exponent for multi-direction flow p:
p = -1.d0 (or any negative value for p) corresponds to an exponent that varies with slope according to: p = 0.5 + 0.6*slope
However, looking at the code in FastScape_ctx.f90:

442     p = pp
443     p_mfd_exp(1:nn) = pp
444     SingleFlowDirection = .false.
445     if (pp.lt.-1.5d0) then
446       SingleFlowDirection = .true.
447       p = 1.d0
448     endif

Here it seems that if the user-set p (pp) is less than -1.5, single flow is assumed and p is set to 1. When SingleFlowDirection is true, then p_mfd_exp is not used (FlowRoutingSingleFlowDirection is called, which does not use p_mfd_exp, instead of FlowRouting, which does use p_mfd_exp which is set to
p_mfd_exp(ij) = 0.5 + 0.6*slope
in FlowRouting.f90.

This makes me think that the description should be
p = -1.5d0 (or any negative value for p smaller than -1.5) corresponds to an exponent of 1, i.e. single direction flow

I get an error building from source on Ubuntu 20.04. Fortran build works without issue, only the python package build fails. I tried the build with --no-build-isolation and got the same error. Attached is a screen dump of the build. The following versions were used.
Python 3.8.10
numpy 1.22.2
wheel 0.37.1
setuptools 60.8.1
scikit-build 0.13.1
cmake 3.22.2
ninja 1.10.2.3

build_screen_dump.txt

As a follow-up on #9, it would be nice to split the following algorithms into separate routines, to allow more flexibility in downstream libraries like fastscape:

• drainage area (currently in StreamPowerLaw)
• closed depression resolving (it has already its own LocalMinima routine, but it would be better to call it in FastScape_Execute_Step instead of inside FlowRouting and FlowRoutingSingleFlowDirection, so that the two latter can be used independently, e.g., from fastscape Python).

This line

Related to #30.

It should be moved in FlowRouting and FlowRoutingSingleFlowDirection, respectively.

Currently all errors are coded using stop. When using the library from Python, this causes the interpreter to end, whereas it would be better to raise an error.

Some related discussion: jameskermode/f90wrap#36

I noticed that the call to the compaction routine is commented out. Is it turned off or is in done somewhere else?

Hello, I've been trying various methods to install your repository (library), but I am encountering an error:

ERROR: Could not find a version that satisfies the requirement fastscapelib_fortran (from versions: none)
ERROR: No matching distribution found for fastscapelib_fortran

I have tried installing through pip, downloading the repository and then installing it, and also through building fastscape.

My system details are:
pip 23.0.1
python 3.7.7
Windows 10 22H2 19045.2604

I am not sure if this information is helpful, but I am hoping you can assist me in resolving this issue.

Currently catchment numbers are computed within the find_stack subroutine, which is called before LocalMinima in flow routing routines.

I think that in most cases we want catchment boundaries to rather correspond to the flow route paths obtained after resolving closed depressions.

@jeanbraun, shouldn't advection affect the vertical elevation of both the bedrock surface and topographic surface with the exact same amount? (or maybe I miss some boundary effects that would cause differences?)

If yes, we might skip solving advection (which is not free) for either h or b and instead apply the delta computed for one surface to the other. We may also get rid of b=min(b,h).

Currently a pyproject.toml (PEP518) is used to specify the tools needed to build the package. This is nice because users don't need to worry about installing these tools first.

However, numpy is needed to both build (using f2py) and import this package at runtime. Build isolation (PEP517) causes issues here. If we add numpy to the build requirements in pyproject.toml, we might have numpy version mismatch at runtime (e.g., when using conda environments) and the package is not usable. If we don't add numpy to the build requirements, we might get a CMake error (numpy not found).

A workaround is to deactivate build isolation, i.e.,

$ pip install . --no-build-isolation

But then we loose the benefits of build requirements in pyproject.toml and users need to install all these tools manually, e.g., using pip or conda:

$ pip install setuptools wheel cmake scikit-build ninja numpy

$ conda install setuptools wheel cmake scikit-build ninja numpy -c conda-forge

Maybe related discussion: pypa/pip#6144

(reported by Rob Moucha)

Currently this maximum number of iterations is hard-coded and cause the solver to stop (SPL), or the whole program to stop (marine):

Users should be able to configure this. We could also allow no limit at all (e.g., a value of -1).

The Union subroutine might cause issues when building the library as a Python module using f2py (see conda-forge/staged-recipes#8461). union is actually a reserved keyword in C/C++ and therefore is handled specifically by f2py.

While this issue might be better solved in f2py, a simple workaround is to rename the (internal) subroutine to avoid conflicts.

@jeanbraun I suspect there is a bug in flexure2D.f90 with the boundary conditions (ibc).

The example flexure_test runs fine:

But when I flip the settings (i.e., setting ibc=0010 and flipping also initial elevation and uplift rates), this is what I get:

The latter doesn't look cyclic on left and right. Erosion is still consistent so I think there might be something wrong in flexure rather than in the model setup.

@benbovy, I have found a memory/assignment error when running the version with split StreamPowerLaw and FlowRouting routines, which I cannot figure out. I propose to rewrite this branch from scratch being more careful on how to split the routine. By creating two versions of the StreamPowerLaw (one for single donor flow and one for multiple donor flow) we should increase clarity and most importantly speed (in the case of single donor flow). I will soon create a new branch with these modifications at which point we should abandon (erase) this branch.

Follow-up on #24, #25.

Now that the Python bindings are available as conda packages.

From Xiaoping Yuan's email:

I see the FastScapeLib Documentation in "https://fastscape-lem.github.io/fastscapelib-fortran/[1]", we shall use O(N) rather than O(n), because n makes readers confused with stream power law exponent n. In the equation, we use capital letter G for the deposition, but in the latter documenting, the little letter g is used. We shall also discuss the equation we put in the document because we have a new version according to our revised JGR paper. The letter p is used to express both slope exponent for multi-direction flow and FastScape_Set_Precip (p) of precipitation rate.

Currently the StreamPowerLaw() routine does many things, actually much more than solving the SPL (i.e., flow routing, closed depression resolving and even applying uplift). This causes various issues:

• hard to read (for example, it took me some time to figure out that uplift is applied inside this routine)
• some potential performance issues (e.g., flow receivers are re-computed in the compute_fluxes routine -- could we avoid that?)
• I have to do hacky things to wrap it with xarray-simlab and make it more modular (see fastscape-lem/fastscape#1)

It would be nice if we could split it into smaller pieces with a clearer scope. In order to still constrain the workflow on the public API, we should also update the FastScape_Execute_Step routine.

cc @jeanbraun