fastscape-lem / fastscapelib-fortran

A Fortran (+ Python bindings) library of efficient algorithms for landscape evolution modeling
https://fastscape-lem.github.io/fastscapelib-fortran/
GNU General Public License v3.0
23 stars 17 forks source link

decrease tolerance and move minima filling code #50

Closed jeanbraun closed 1 year ago

jeanbraun commented 1 year ago

@benbovy Could you please accept this pull request and rebuild the conda version/depot of Fastscape so that others (and especially Amanda Wild) can use this new version. Note that Amanda is waiting for this new version to finish a couple of figures on sensitivity analysis for her paper. There are two issues resolved with this pull request: 1) solves a problem of mass conservation noted by Sebastian Wolf when several minima occur along a flow line; for this, as suggested by Sebastian, the code performing the sediment fill of the local minima has been moved to the end of the StreamPowerLaw.f90 routine (both for the single and multiple flow algorithms) 2) solves an issue when performing very small time steps that required the tolerance to be reduced by a couple of orders of magnitude; note that this is not a proper fix but I am currently working on an updated version of the multi receiver algorithm with Guillaume Cordonnier in which we will pay more attention to this convergence criterion/tolerance Finally, note that all changes were performed in the StreamPowerLaw.f90 file and fully tested using a Fortran test as I could not recompile the python bidding on my computer (M1 running Venture 13.0.1). Thanks in advance.

benbovy commented 1 year ago

Thanks @jeanbraun, I'll do a new release and build new conda-forge packages (python bindings).

I just want to check with you first if the new updated value of tolerance won't have a major impact on users of fastscapelib-fortran or fastscape (python). Since the tolerance is hard-coded it is impossible to change its value, and I also see that you removed the max. number of iterations. There's a risk of having infinite loop if it doesn't converge, do you think the tolerance value is safe enough here?

cc @rbrtmch who reported on Gitter some issues with this using specific values of precipitation rate and g.

Ideally we should expose both the tolerance and nb. of iterations as parameters. See #44.

jeanbraun commented 1 year ago

Hello @benbovy. I fully agree with you that we should avoid breaking changes for all users. I note, however, that the condition on the maximum number of iterations (100) is still there but it does not produce the warning message that was cluttering the output of some runs. I would be happy to reinstate it if necessary. But the algorithm should never get into an infinite loop. I also agree that it would be ideal to expose both the tolerance and the maximum number of iterations permitted. I am happy to do this. Let me know.

benbovy commented 1 year ago

I note, however, that the condition on the maximum number of iterations (100) is still there

Ah yes you're right, sorry I've seen the diffs too quickly.

It think it's fine to remove the message (not sure it was displayed with the Python bindings anyway).

If we assume that >100 iterations would occur on very rare occasions with the new tolerance value, I think we can go ahead with the release if that's urgent. I can still merge Sebastian's pull-request and we could do the same for the max. number of iterations if that's necessary.