Fix updating of the topographic elevations

news/67.bugfix

@@ -0,0 +1,3 @@
+
+Fixed a bug where the topographic elevations were not being updated correctly after the
+components were time stepped.

news/67.feature

@@ -0,0 +1,4 @@
+
+Added a :meth:`~sequence.SequenceModelGrid.get_profile` method to :class:`~sequence.SequenceModelGrid` that is a convenience
+function to return a field just along the middle row of nodes (i.e. along the
+profile).

news/67.feature.1

@@ -0,0 +1,6 @@
+
+Changed all components to only update the current amount of subsidence and
+deposited sediment but not to update secondary fields like the topographic
+elevation. This removes the dependence of these secondary fields on
+component ordering and ensures these secondary fields are updated once, and only
+once, with each time step.

news/67.misc

@@ -0,0 +1,2 @@
+
+Added tests for the compaction, flexure, and subsidence components.

sequence/_grid.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 from landlab import RasterModelGrid
+from numpy.typing import NDArray
 
 if sys.version_info >= (3, 11):
     import tomllib
@@ -41,6 +42,21 @@ def __init__(self, n_cols: int, spacing: float = 100.0):
         self.at_node["sediment_deposit__thickness"] = self.zeros(at="node")
         self.at_grid["sea_level__elevation"] = 0.0
 
+    def get_profile(self, name: str) -> NDArray:
+        """Return the values of a field along the grid's profile.
+
+        Parameters
+        ----------
+        name : str
+            The name of an at-node field.
+
+        Returns
+        -------
+        values : ndarray
+            The values of the field located at the middle row of nodes.
+        """
+        return self.at_node[name].reshape(self.shape)[1]
+
     @classmethod
     def from_toml(cls, filepath: os.PathLike[str]) -> "SequenceModelGrid":
         """Load a :class:`~SequenceModelGrid` from a *toml*-formatted file.

sequence/cli.py

@@ -231,8 +231,12 @@ def run(ctx: Any, with_citations: bool, dry_run: bool) -> None:
     run_dir = pathlib.Path.cwd()
 
     times, names = _find_config_files(".")
+    if len(times) == 0:
+        err("[ERROR] unable to find a configuration file.")
+        raise click.Abort()
+
     if not silent:
-        out(f"config_files = [{', '.join(repr(name) for name in names)}]")
+        out(f"[INFO] config_files: {', '.join(repr(name) for name in names)}")
     params = TimeVaryingConfig.from_files(names, times=times)
 
     model_params = params.as_dict()
@@ -250,20 +254,16 @@ def run(ctx: Any, with_citations: bool, dry_run: bool) -> None:
     )
 
     if verbose or not silent:
-        out("enabled = [")
         for name in model.components:
-            out(f"  {name!r},")
-        out("]")
-        out("disabled = [")
+            out(f"[INFO] ✅ enabled: {name}")
         for name in set(SequenceModel.ALL_PROCESSES) - set(model.components):
-            out(f"  {name!r},")
-        out("]")
+            out(f"[INFO] ❌ disabled {name}")
 
     if not silent and verbose:
         out(params.dump())
 
     if not silent and len(processes) == 0:
-        out("⚠️  ALL PROCESSES HAVE BEEN DISABLED! ⚠️")
+        out("[WARN] ⚠️  ALL PROCESSES HAVE BEEN DISABLED! ⚠️")
 
     if not silent and with_citations:
         from landlab.core.model_component import registry

sequence/fluvial.py

@@ -34,6 +34,14 @@ class Fluvial(Component):
             "mapping": "grid",
             "doc": "Position of sea level",
         },
+        "sediment_deposit__thickness": {
+            "dtype": "float",
+            "intent": "inout",
+            "optional": False,
+            "units": "m",
+            "mapping": "node",
+            "doc": "Thickness of deposition or erosion",
+        },
         "delta_sediment_sand__volume_fraction": {
             "dtype": "float",
             "intent": "out",
@@ -135,17 +143,17 @@ def run_one_step(self, dt: float) -> None:
         channel_width = sand_vol / qs
 
         x = self.grid.x_of_node.reshape(self.grid.shape)[1]
-        z = self.grid.at_node["topographic__elevation"].reshape(self.grid.shape)[1]
+        elevation = self.grid.get_profile("topographic__elevation")
 
         shore = find_shoreline(
-            x, z, sea_level=self.grid.at_grid["sea_level__elevation"]
+            x, elevation, sea_level=self.grid.at_grid["sea_level__elevation"]
         )
 
         land = x < shore
         water = x >= shore
         # land = self.grid.x_of_node[self.grid.node_at_cell] < shore
         # slope = np.gradient(z[1, land]) / self.grid.dx
-        slope = np.gradient(z) / self.grid.dx
+        slope = np.gradient(elevation) / self.grid.dx
         # slp = np.zeros(1)
         # slp[0] = self.plain_slope
         # slope = np.concatenate((slp,slop))
@@ -155,7 +163,7 @@ def run_one_step(self, dt: float) -> None:
         # )
 
         # channel_depth = np.zeros(self.grid.number_of_cells)
-        channel_depth = np.zeros_like(z)  # use upsteam slope
+        channel_depth = np.zeros_like(elevation)  # use upsteam slope
         # channel_depth[0] = (
         #    (self.sand_density - 1000.) / 1000. * self.sand_grain / self.plain_slope
         channel_depth[land] = (
@@ -170,28 +178,13 @@ def run_one_step(self, dt: float) -> None:
         # width_cb = channel_width/epsilon
 
         # Original: r_cb = (model.new_height[i]-model.height[i]+model.d_sl);
-        r_cb = (
-            self.grid.at_node["sediment_deposit__thickness"]
-            .reshape(self.grid.shape)[1]
-            .copy()
-            * epsilon
-        )
-        dz = (
-            self.grid.at_node["bedrock_surface__increment_of_elevation"]
-            .reshape(self.grid.shape)[1]
-            .copy()
-        )
+        r_cb = self.grid.get_profile("sediment_deposit__thickness") * epsilon
+        dz = self.grid.get_profile("bedrock_surface__increment_of_elevation").copy()
         # original: r_b = model.thickness[i];
-        r_b = (
-            self.grid.at_node["sediment_deposit__thickness"]
-            .reshape(self.grid.shape)[1]
-            .copy()
-        )
+        r_b = self.grid.get_profile("sediment_deposit__thickness").copy()
         # r_fp = np.zeros(self.grid.shape[1])
-        r_fp = np.zeros_like(z)
-        percent_sand = self.grid.at_node[
-            "delta_sediment_sand__volume_fraction"
-        ].reshape(self.grid.shape)[1]
+        r_fp = np.zeros_like(elevation)
+        percent_sand = self.grid.get_profile("delta_sediment_sand__volume_fraction")
         percent_sand.fill(0.0)
 
         for i in np.where(land)[0]:
@@ -295,25 +288,13 @@ def run_one_step(self, dt: float) -> None:
                 if add_mud[i] < 0.0:
                     add_mud[i] = 0.0
 
-        thickness = (
-            self.grid.at_node["sediment_deposit__thickness"]
-            .reshape(self.grid.shape)[1]
-            .copy()
-        )
-        thickness[water] += add_mud[water]
+        sediment_thickness = self.grid.get_profile("sediment_deposit__thickness")
+        sediment_thickness[water] += add_mud[water]
 
         np.divide(
-            thickness[water] - add_mud[water],
-            thickness[water],
-            where=thickness[water] > 0.0,
+            sediment_thickness[water] - add_mud[water],
+            sediment_thickness[water],
+            where=sediment_thickness[water] > 0.0,
             out=percent_sand[water],
         )
-
         np.clip(percent_sand[water], 0.0, 1.0, out=percent_sand[water])
-
-        plus_mud = np.zeros(self.grid.shape)
-        plus_mud[1] = add_mud
-        sdt = self.grid.at_node["sediment_deposit__thickness"].reshape(plus_mud.shape)
-        sdt[1][:] += plus_mud[1]
-        te = self.grid.at_node["topographic__elevation"].reshape(plus_mud.shape)
-        te[1][:] -= plus_mud[1]