Merge pull request #221 from amoodie/joss_salterg_review

Salter JOSS review reply - part2

docs/source/guides/10min.rst

@@ -13,7 +13,7 @@ This simple guide will show you the absolute basics of getting a `pyDeltaRCM` mo
 A default model
 ---------------
 
-You can get a model running with five simple lines of code.
+You can get a model running with five simple lines of Python code. Note that you can run *pyDeltaRCM* in either a standalone script or part of an interactive session.
 First, we instantiate the main :obj:`~pyDeltaRCM.DeltaModel` model object.
 
 .. code:: python

docs/source/guides/advanced_configuration_guide.inc

@@ -34,7 +34,7 @@ Each folder will contain a copy of the configuration file used for that job; for
 .. code:: yaml
 
     out_dir: 'out_dir/job_000'
-    dx: 2.0
+    dx: 40.0
     h0: 1.0
     f_bedload: 0.5
 

docs/source/guides/developer_guide.rst

@@ -214,7 +214,20 @@ From these results, we can see that the values returned from the built-in unifor
 Model development
 -----------------
 
-.. todo::
+Slicing and neighbors 
+~~~~~~~~~~~~~~~~~~~~~
 
-    add some notes about slicing arrays, and how we pad them with a custom pad operation when needed. Look at shared tools docs for starting point.
+Slicing an array to find the array values of neighbors is a common operation in the model.
+The preferred way to slice is by 1) padding the array with :func:`~pyDeltaRCM.shared_tools.custom_pad`, and 2) looping through rows and columns to directly index. This approach makes for readable and reasonably fast code; for example, to find any cells that are higher than all neighbors:
 
+.. code::
+    
+    pad_eta = shared_tools.custom_pad(self.eta)
+    for i in range(self.L):
+        for j in range(self.W):
+            eta_nbrs = pad_eta[i - 1 + 1:i + 2 + 1, j - 1 + 1:j + 2 + 1]
+            eta_nbrs[1, 1] = -np.inf
+            
+            np.all(self.eta[i, j] > eta_nbrs)
+
+There are also several model attributes that may be helpful in development; we suggest using these builtins rather than creating your own whenever possible (see :meth:`~pyDeltaRCM.init_tools.init_tools.set_constants` and the model source code).

docs/source/info/morphodynamics.rst

@@ -18,18 +18,34 @@ In this documentation, we focus on the details of *model implementation*, rather
 Sediment Transport
 ==================
 
-.. note::
-   Incomplete.
+Sediment transport in the model is computed according to an excess stress approach.
+Conceptually, sand is routed as bed-material load, and mud is routed as fully suspended load.
+
+For sand parcels, the *transport capacity* is determined by the scaling between sediment flux and flow velocity, and takes the form of the Meyer-Peter and Müller (1948) [3]_ formula:
+
+.. math::
+
+      q_{s\_cap} = q_{s0} \frac{u^\beta_{loc}}{u^\beta_0},
+
+where :math:`u_{loc}` is the depth averaged flow velocity in the cell, :math:`beta` is an exponent set to 3 by default (:obj:`~pyDeltaRCM.model.DeltaModel.beta`), and :math:`q_{s0}` is the unit-width upstream sand flux input at the inlet channel.
+At each step of the model domain, sand is either eroded or deposited to the bed depending on the local flow velocity :math:`u_{loc}` and local sediment transport :math:`q_{s\_loc}`. 
+Sand is deposited where local transport (i.e., the sand put into that cell from upstream) is greater than the cell transport capacity :math:`q_{s\_loc} > q_{s\_cap}`.
+Sand is eroded from the bed when the local velocity is greater than the threshold erosion velocity (:obj:`~pyDeltaRCM.model.DeltaModel.coeff_U_ero_sand`) **and** the local transport is less than the local transport capacity.
+
+
+Mud parcels do not have any local capacity (i.e., fully suspended washload transport). 
+At each parcel step, mud is either eroded or deposited (or neither), depending on the relative value of local flow velocity :math:`u_{loc}` and the threshold erosion and deposition values (:obj:`~pyDeltaRCM.model.DeltaModel.coeff_U_ero_mud` and :obj:`~pyDeltaRCM.model.DeltaModel.coeff_U_dep_mud`).
+
+.. note:: 
+
+      A complete conceptual description of sediment erosion and deposition routing rules can be found in the original DeltaRCM reference Liang et al., 2015 [1]_.
 
 
 .. _sediment-routing-weighting:
 
 Sediment routing weighting
 --------------------------
 
-.. note::
-   Incomplete.
-
 Sediment routing probability for a given cell :math:`j` to neighbor cell :math:`i` is computed according to:
 
 .. math::
@@ -65,7 +81,7 @@ where :math:`\Delta V` is the volume of sediment to be eroded or deposited from
 
 .. note::
 
-    Total sediment mass is preserved, but individual categories of sand and mud are not. I.e., it is assumed that there is an infinite supply of sand and/or mud to erode and entrain at any location in the model domain.
+    Total sediment mass is preserved during erosion, but individual categories of sand and mud are not. I.e., it is assumed that there is an infinite supply of sand and/or mud to erode and entrain at any location in the model domain.
 
 Following a change in the bed elevation, the local flow depth is updated and then local flow velocity is updated according to fluid mass conservation (i.e., ``uw = qw / h``; :obj:`~sed_tools.BaseRouter._update_fields`; [1]_).
 
@@ -158,8 +174,8 @@ The impact of topographic diffusion is minor compared to the bed elevation chang
 Notes for modeling best practices
 =================================
 
-* Stop simulations before the delta reaches the edge of the computational domain. If a distributary channel reaches the domain edge, this channel is likely to become locked in place, and will convey sediment outside the computational domain, thus violating any statements of mass conservation. Generally, simulations that reach the edge of the domain should be discarded. 
-* Stop simulations before the delta reaches a grade condition, where the topset slope is equal to background slope `S0`. This is really only an issue for large domains run for long duration.
+* Stop simulations before the delta reaches the edge of the computational domain. Delta channel dynamics are changed when a distributary channel reaches the domain edge, because sediment is conveyed to outside the computational domain where it no longer feeds back on channel development. Channels become "locked" in place in this scenario [2]_, because the domain edge is an infinite sediment sink, and therefore rendering invalid any assumptions about stationarity of delta dynamics and/or stratigraphy. Moreover, the downstream water surface boundary condition (`H_SL`) will be violated if a channel reaches the domain edge. Generally, simulations that reach the edge of the domain should be discarded.
+* Stop simulations before the delta reaches a condition where the topset slope is equal to the background slope parameter (`S0`). When the background slope is reached, the transport capacity of sediment through the delta is diminished such that channels "clog" up and trigger model instabilities. This is really only an issue for large domains run for long duration.
 * Use a sufficient number of water and sediment parcels (> 2000). Too few parcels will result in a rough water surface, irregular sediment deposition, and a rough bed elevation.
 
 
@@ -170,3 +186,9 @@ References
 .. [1] A reduced-complexity model for river delta formation – Part 1: Modeling
        deltas with channel dynamics, M. Liang, V. R. Voller, and C. Paola, Earth
        Surf. Dynam., 3, 67–86, 2015. https://doi.org/10.5194/esurf-3-67-2015
+
+.. [2] Liang, M., Kim, W., and Passalacqua, P. (2016), How much subsidence is
+       enough to change the morphology of river deltas?, Geophysical Research Letters, 43, 10,266--10,276, doi:10.1002/2016GL070519.
+
+.. [3] Meyer-Peter, E. and Müller, R.: Formulas for bed-load transport, in: 
+       Proceedings of the 2nd Meeting of IAHSR, Stockholm, Sweden, 39–64, 1948.

docs/source/meta/contributing.rst

@@ -2,17 +2,28 @@
 Contributing
 ************
 
-This project follows the :doc:`Code of Conduct </meta/conduct>`. 
+To contribute to this project, you must follow our :doc:`Code of Conduct </meta/conduct>` at all times.
 If you are not familiar with our code of conduct policy, take a minute to read the policy before starting with your first contribution.
 
+
 How to contribute
 -----------------
 
-If you are interested in contributing, please submit a pull request or get in touch with the development team via the `Github issue tracker <https://github.com/DeltaRCM/pyDeltaRCM/issues>`_.
-
 We welcome contributions of many types, including but not limited to:
 
     * bug fixes
     * improvements to documentation
     * new features
     * additional examples for using and developing *pyDeltaRCM*
+
+If you are interested in contributing, please submit a pull request or get in touch with the development team via the `Github issue tracker <https://github.com/DeltaRCM/pyDeltaRCM/issues>`_.
+
+
+Issues and questions
+--------------------
+
+If you have identified a bug within the code, but aren't sure how to (or don't want to) fix it, we invite you to open an issue on our `Github issue tracker <https://github.com/DeltaRCM/pyDeltaRCM/issues>`_.
+Please explain the problem as fully as possible, including any system information that may be relevant.
+
+You can also ask any questions about the software or ask for help on the `Github issue tracker <https://github.com/DeltaRCM/pyDeltaRCM/issues>`_.
+We will try our best to help you how we can!

docs/source/meta/installing.rst

@@ -26,7 +26,7 @@ With Anaconda on Linux:
     $ conda activate deltarcm
 
 For more informtaion, see `this guide <https://conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html#activating-an-environment>`_ for help on creating and activating a virtual environment with Anaconda on other platforms.
-See `this guide <https://packaging.python.org/guides/installing-using-pip-and-virtual-environments/#creating-a-virtual-environment>`_ for creating virtual environments with `venv` if you do not use Anaconda.
+See `this helpful guide <https://packaging.python.org/guides/installing-using-pip-and-virtual-environments/#creating-a-virtual-environment>`_ for creating virtual environments with `venv` if you do not use Anaconda.
 
 
 User installation

docs/source/reference/init_tools/index.rst

@@ -4,9 +4,33 @@
 init_tools
 **********
 
-.. todo::
+.. currentmodule:: pyDeltaRCM.init_tools
+
+
+The model initialization is managed by :obj:`~pyDeltaRCM.model.DeltaModel.__init__`, but the actual initialization is mostly handled by routines in `init_tools`.
+The major steps of initialization are:
+
+.. autosummary::
+
+    init_tools.import_files
+    init_tools.init_logger
+    init_tools.process_input_to_model
+    init_tools.determine_random_seed
+    init_tools.create_other_variables
+    init_tools.create_domain
+    init_tools.init_sediment_routers
+    init_tools.init_subsidence
+
+and then depending on the checkpointing configuration, the following methods may be called:
+
+.. autosummary::
+
+    init_tools.load_checkpoint
+    init_tools.init_output_file
+    pyDeltaRCM.interation_tools.output_data
+    pyDeltaRCM.interation_tools.output_checkpoint
+    pyDeltaRCM.interation_tools.log_model_time
 
-    Add paragraph description of the module. What stages are defined here generally? Make a table with the main ones like in water tools?
 
 Public API methods attached to model
 ------------------------------------

pyDeltaRCM/init_tools.py

@@ -283,6 +283,16 @@ def set_constants(self):
         Configure constants, including coordinates and distances, as well as
         environmental constants (gravity), and kernels for smoothing
         topography.
+
+        Some of the constants defined herein:
+            * `self.g`, gravitational acceleration
+            * `self.distances`, distance from cell `i,j` to neighbors (and self)  # noqa: E501
+            * `self.iwalk`, step distance cross domain to cell in indexed direction
+            * `self.jwalk`, step distance down domain to cell in indexed direction
+            * `self.ravel_walk`, flattened index distance to cell in indexed direction
+
+        Each of these attributes also has a `self.xxxxx_flat` sibling
+        attribute, which is simply the flattened version of that attribute.
         """
         _msg = 'Setting model constants'
         self.log_info(_msg, verbosity=1)

pyDeltaRCM/shared_tools.py

@@ -150,8 +150,26 @@ def get_random_uniform(limit):
 def get_start_indices(inlet, inlet_weights, num_starts):
     """Get start indices.
 
-    .. todo:: Description needed. Where is it used?
+    Reutrn a randomly generated list of starting points for parcel routing.
+    These starting points are selected from the `inlet` array.
 
+    Parameters
+    ----------
+    inlet : ndarray
+        Array of inlet cells.
+
+    inlet_weights : ndarray
+        Array of weights to select items from :obj:`inlet`. Should have same
+        dimensions.
+
+    num_starts : int
+        Number of starting points to generate.
+
+    Returns
+    -------
+    start_indices : ndarray
+        :obj:`num_starts` starting points, generated from :obj:`inlet`
+         according to weights in :obj:`inlet_weights`. 
     """
     norm_weights = inlet_weights / np.sum(inlet_weights)
     idxs = []
@@ -162,14 +180,14 @@ def get_start_indices(inlet, inlet_weights, num_starts):
 
 
 @njit
-def get_steps(new_cells, iwalk, jwalk):
-    """Find the values giving the next step.
-
-    .. todo:: Description needed. Where is it used?
+def get_steps(new_direction, iwalk, jwalk):
+    """Find the values given the next step.
 
+    Get the steps for updating discharge and velocity arrays based on the
+    direction of each step.
     """
-    istep = iwalk[new_cells]
-    jstep = jwalk[new_cells]
+    istep = iwalk[new_direction]
+    jstep = jwalk[new_direction]
     dist = np.sqrt(istep * istep + jstep * jstep)
     astep = (dist != 0)