Merge pull request #172 from amoodie/docs_sediment

Docs water and sediment routing (partial)

docs/source/info/hydrodynamics.rst

@@ -11,15 +11,53 @@ In this documentation, we focus on the details of *model implementation*, rather
 Routing individual water parcels
 ================================
 
-.. note::
+Probabilities for water parcel routing *to all neighbors and to self* for *each cell* are computed *once* at the beginning of the water routing routine (:obj:`~water_tools.get_water_weight_array` called from :obj:`~water_tools.run_water_iteration`).
 
-   Incomplete.
+Water routing probability for a given cell :math:`j` to neighbor cell :math:`i` is computed according to:
 
+.. math::
+
+    w_i = \frac{\frac{1}{R_i} \max(0, \mathbf{F}\cdot\mathbf{d_i})}{\Delta i},
+
+where :math:`\mathbf{F}` is the local routing direction and :math:`\mathbf{d_i}` is a unit vector pointing to neighbor :math:`i` from cell :math:`j`, and :math:`\Delta_i` is the cellular distance to neighbor :math:`i` (:math:`1` for cells in main compass directions and :math:`\sqrt{2}` for corner cells.
+:math:`R_i` is a flow resistance estimated as an inverse function of local water depth (:math:`h_i`):
+
+.. math::
+
+    R_i = \frac{1}{{h_i}^\theta}
+
+The exponent :math:`\theta` takes a value of unity by default for water routing (:attr:`~pyDeltaRCM.DeltaModel.theta_water`, :doc:`/reference/model/yaml_defaults`), leading to routing weight for neighbor cell :math:`i`:
+
+.. math::
+
+    w_i = \frac{h_i \max(0, \mathbf{F}\cdot\mathbf{d_i})}{\Delta i},
+
+These weights above are calculated only for wet neighbor cells; all dry neighbor cells take a weight value of 0 (:obj:`~water_tools._get_weight_at_cell_water`).
+Finally, probability for routing from cell :math:`j` to cell :math:`i` is calculated as:
+
+.. math::
+
+    p_i = \frac{w_i}{\sum^8_{nb=1} w_{nb}}, i=1, 2, \ldots, 8
+
+Weights are accumulated for 8 neighbors and a probability of 0 is assigned to moving from cell :math:`j` to cell :math:`j` (i.e., no movement).
+These 9 probabilities are organized into an array ``self.water_weights`` with shape (:obj:`L`, :obj:`W`, 9)`. 
 
 The following figure shows several examples of locations within the model domain, and the corresponding water routing weights determined for that location.
 
 .. plot:: water_tools/water_weights_examples.py
 
+Because probabilities are computed for all locations once at the beginning of water iteration, all water parcels can be routed *in parallel* step-by-step in :obj:`~water_tools.run_water_iteration`.
+During iteration, the direction of the random walk is chosen for each parcel via :obj:`_choose_next_directions`, which internally uses :func:`~pyDeltaRCM.shared_tools.random_pick` for randomness.
+For example, see the random walks of several parcels below:
+
+.. plot::  water_tools/run_water_iteration.py
+
+
+.. todo:: add sentence or two above about check_for_loops.
+
+Water routing completes when all water parcels have either 1) reached the model domain boundary, 2) taken a number of steps exceeding :attr:`~pyDeltaRCM.model.DeltaModel.stepmax`, or 3) been removed from further routing via the :obj:`_check_for_loops` function.
+
+
 
 Combining parcels into free surface
 ===================================
@@ -45,16 +83,18 @@ The updated water surface is combined with the previous timestep's water surface
 
 With a new free surface computed, a few final operations prepare the surface for boundary condition updates and eventually being passed to the sediment routing operations.
 A non-linear smoothing operation is applied to the free surface, whereby wet cells are iteratively averaged with neighboring wet cells to yield an overall smoother surface.
-The smoothing is handled by :func:`_smooth_free_surface` and depends on the number of iterations (:attr:`~pyDeltaRCM.DeltaModel.Nsmooth`) and a weighting coefficient (:attr:`~pyDeltaRCM.DeltaModel.Csmooth`).
+The smoothing is handled by :func:`_smooth_free_surface` and depends on the number of iterations (:attr:`~pyDeltaRCM.model.DeltaModel.Nsmooth`) and a weighting coefficient (:attr:`~pyDeltaRCM.model.DeltaModel.Csmooth`).
 
 .. plot:: water_tools/_smooth_free_surface.py
 
 
-.. todo:: add component describing the smoothing.
-
-.. todo:: add component describing the flooding correction.
+Finally, a :meth:`~water_tools.flooding_correction` is applied to the domain.
+In this correction, all "dry" cells (a cell where the flow depth is less than the `dry_depth`) are checked for any neighboring cells where the water surface elevation (`stage`) is higher than the bed elevation of the dry cell.
+If this condition is met for a given dry cell, the dry cell is flooded: the stage of the dry cell is set to the maximum stage of neighboring cells.
 
+.. plot:: water_tools/flooding_correction.py
 
+Similar to :func:`~_smooth_free_surface` described above, :meth:`~water_tools.flooding_correction` acts to remove roughness in the water surface and contributes to :ref:`model stability <model-stability>`.
 
 
 Finalizing and boundary conditions to sediment routing

docs/source/info/morphodynamics.rst

@@ -12,6 +12,8 @@ See [1]_ for a complete description of morphodynamic assumptions in the DeltaRCM
 In this documentation, we focus on the details of *model implementation*, rather than *model design*.
 
 
+.. _sediment-transport:
+
 ==================
 Sediment Transport
 ==================
@@ -20,24 +22,69 @@ Sediment Transport
    Incomplete.
 
 
+.. _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::
+
+    w_i = \frac{\frac{1}{R_i} \max(0, \mathbf{F}\cdot\mathbf{d_i})}{\Delta i},
+
+where :math:`\mathbf{F}` is the local routing direction and :math:`\mathbf{d_i}` is a unit vector pointing to neighbor :math:`i` from cell :math:`j`, and :math:`\Delta_i` is the cellular distance to neighbor :math:`i` (:math:`1` for cells in main compass directions and :math:`\sqrt{2}` for corner cells.
+:math:`R_i` is a resistance estimated as an inverse function of local water depth (:math:`h_i`):
+
+.. math::
+
+    R_i = \frac{1}{{h_i}^\theta}.
+
+Here, :math:`\theta` takes the value of :obj:`~pyDeltaRCM.model.DeltaModel.coeff_theta_sand` for sand routing probabilities, and :obj:`~pyDeltaRCM.model.DeltaModel.coeff_theta_mud` for mud routing.
+
+
+.. plot:: sed_tools/sediment_weights_examples.py
+
+
 ============================
 Changes in the bed elevation
 ============================
 
-Change in the channel bed is the result of deposition or erosion by each sediment parcel (:obj:`_update_fields`), dictated by sediment mass conservation (i.e., Exner equation) and is equal to:
+Along the walk of a sediment parcel, the sediment parcel volume is modulated on each step, according to the sediment transport rules described above in :ref:`sediment-transport`.
+As the volume of the sediment parcel changes, the channel bed elevation at the current parcel location is updated to reflect this volume change (:obj:`~sed_tools.BaseRouter._update_fields`), i.e., the bed is eroded or sediment is deposited on the bed.
+The vertical change in the channel bed is dictated by sediment mass conservation (i.e., Exner equation) and is equal to:
 
 .. math::
 
     \Delta \eta = \Delta V / dx^2
 
+where :math:`\Delta V` is the volume of sediment to be eroded or deposited from the bed at a given cell along the parcel walk.
 
 .. note::
 
-    Total sediment mass is preserved, but individual categories of sand and mud are not. I.e., it is assumed that there is infinite sand and/or mud to erode at any location where erosion is occurring.
+    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.
 
-.. todo::
+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]_).
 
-   Incomplete.
+Sediment parcels are routed through the model domain step-by-step and in serial, such that changes in the bed elevation caused by one sediment parcel will affect the weighted random walks of all subsequent sediment parcels (:ref:`sediment-routing-weighting`), due to the updated flow field.
+
+Sediment parcel routing is handled by first routing all sand parcels, applying a topographic diffusion (see below and :meth:`~sed_tools.topo_diffusion`), and then routing all mud parcels.
+The impact of routing *all* sand and mud parcels on bed elevation is shown in the table below.
+
+.. _sand-mud-route-comparison:
+
+.. table::
+
+    +-------------------------------------------+-----------------------------------------------+----------------------------------------------+
+    | initial bed                               | :meth:`~sed_tools.route_all_sand_parcels`     | :meth:`~sed_tools.route_all_mud_parcels`     |
+    +===========================================+===============================================+==============================================+
+    | .. plot:: sed_tools/_initial_bed_state.py | .. plot:: sed_tools/route_all_sand_parcels.py | .. plot:: sed_tools/route_all_mud_parcels.py |
+    +-------------------------------------------+-----------------------------------------------+----------------------------------------------+
+
+.. _model-stability:
 
 ===============
 Model Stability
@@ -48,6 +95,25 @@ Model stability depends on...
 .. note::
    Incomplete.
 
+.. topographic-diffusion:
+
+Topographic diffusion
+---------------------
+
+Abrupt change in bed elevation (i.e., steep local bed slope) may lead to numerical instability. 
+To prevent this, a topographic diffusion is applied immediately following the routing of all sand parcels in the model sequence.
+
+.. hint::
+
+    Topographic diffusion is applied between routing sand parcels and routing mud parcels.
+
+In implementation, topographic smoothing convolves topography with `3x3` cell kernels configured to a diffusive behavior.
+The diffusion is repeated over the entire model domain :obj:`~pyDeltaRCM.DeltaModel.N_crossdiff` times.
+In the following example, :obj:`~pyDeltaRCM.DeltaModel.N_crossdiff` takes the :doc:`default value </reference/model/yaml_defaults>`.
+
+.. plot:: sed_tools/topo_diffusion.py
+
+The impact of topographic diffusion is minor compared to the bed elevation change driven by parcel erosion or deposition (:ref:`sand and mud routing effects <sand-mud-route-comparison>`).
 
 .. _reference-volume:
 

docs/source/pyplots/sed_tools/_initial_bed_state.py

@@ -0,0 +1,30 @@
+import warnings
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import pyDeltaRCM
+
+# filter out the warning raised about no netcdf being found
+warnings.filterwarnings("ignore", category=UserWarning)
+
+
+# init delta model
+with pyDeltaRCM.shared_tools._docs_temp_directory() as output_dir:
+    delta = pyDeltaRCM.DeltaModel(
+        out_dir=output_dir,
+        resume_checkpoint='../../_resources/checkpoint')
+
+_shp = delta.eta.shape
+
+# set up axis
+fig, ax = plt.subplots()
+
+# fill in axis2
+pyDeltaRCM.debug_tools.plot_domain(
+    delta.eta, ax=ax, grid=False, cmap='cividis')
+ax.set_title('intial bed elevation (m)')
+
+
+plt.tight_layout()
+plt.show()

docs/source/pyplots/sed_tools/route_all_mud_parcels.py

@@ -0,0 +1,39 @@
+import warnings
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import pyDeltaRCM
+
+# filter out the warning raised about no netcdf being found
+warnings.filterwarnings("ignore", category=UserWarning)
+
+
+# init delta model
+with pyDeltaRCM.shared_tools._docs_temp_directory() as output_dir:
+    delta = pyDeltaRCM.DeltaModel(
+        out_dir=output_dir,
+        resume_checkpoint='../../_resources/checkpoint')
+
+_shp = delta.eta.shape
+
+_eta_before = np.copy(delta.eta)
+
+delta.route_all_mud_parcels()
+
+_eta_after = np.copy(delta.eta)
+
+
+# set up axis
+fig, ax = plt.subplots()
+
+# fill in axis2
+_diff = _eta_after - _eta_before
+pyDeltaRCM.debug_tools.plot_domain(
+    _diff, ax=ax, grid=False, cmap='RdBu',
+    vmin=-0.1, vmax=0.1)
+ax.set_title('bed elevation change (m) \n due to mud parcels')
+
+
+plt.tight_layout()
+plt.show()

docs/source/pyplots/sed_tools/route_all_sand_parcels.py

@@ -0,0 +1,39 @@
+import warnings
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import pyDeltaRCM
+
+# filter out the warning raised about no netcdf being found
+warnings.filterwarnings("ignore", category=UserWarning)
+
+
+# init delta model
+with pyDeltaRCM.shared_tools._docs_temp_directory() as output_dir:
+    delta = pyDeltaRCM.DeltaModel(
+        out_dir=output_dir,
+        resume_checkpoint='../../_resources/checkpoint')
+
+_shp = delta.eta.shape
+
+_eta_before = np.copy(delta.eta)
+
+delta.route_all_sand_parcels()
+
+_eta_after = np.copy(delta.eta)
+
+
+# set up axis
+fig, ax = plt.subplots()
+
+# fill in axis2
+_diff = _eta_after - _eta_before
+pyDeltaRCM.debug_tools.plot_domain(
+    _diff, ax=ax, grid=False, cmap='RdBu',
+    vmin=-0.1, vmax=0.1)
+ax.set_title('bed elevation change (m) \n due to sand parcels')
+
+
+plt.tight_layout()
+plt.show()