Merge branch 'master' into bc/space_return_eff_ero_depo

news/1967.bugfix

@@ -0,0 +1,2 @@
+
+Update ConcentrationTrackerForDiffusion run methods to to solve two situations that cause the mass balance to fail.

src/landlab/components/concentration_tracker/concentration_tracker_for_diffusion.py

@@ -30,12 +30,14 @@ class ConcentrationTrackerForDiffusion(Component):
     .. note::
 
         This component requires a soil flux field calculated by a hillslope diffusion
-        component and must be run after every diffusion step. Currently, this component
-        WILL ONLY WORK IF COUPLED with the :class:`~.DepthDependentDiffuser` or the
-        :class:`~.DepthDependentTaylorDiffuser` (without the dynamic timestep option).
+        component. This component is implemented by applying a :meth:`start_tracking`
+        method immediately before every diffusion step and a :meth:`stop_tracking`
+        method immediately after every diffusion step. These methods are applied instead
+        of the typical :meth:`run_one_step` method. See the docstring examples below.
 
-        In-situ production and decay of the material property are handled by
-        the ConcentrationTrackerProductionDecay component.
+        Currently, this component will only work if coupled with the
+        :class:`~.DepthDependentDiffuser` or the :class:`~.DepthDependentTaylorDiffuser`
+        (without the dynamic timestep option).
 
     Examples
     --------
@@ -75,8 +77,9 @@ class ConcentrationTrackerForDiffusion(Component):
 
     >>> ddd = DepthDependentDiffuser(mg)
     >>> ct = ConcentrationTrackerForDiffusion(mg)
+    >>> ct.start_tracking()
     >>> ddd.run_one_step(1.0)
-    >>> ct.run_one_step(1.0)
+    >>> ct.stop_tracking(1.0)
 
     >>> mg.at_node["topographic__elevation"].reshape(mg.shape)
     array([[ 0. ,  4.        ,  8.        , 12.        , 16. ],
@@ -91,8 +94,11 @@ class ConcentrationTrackerForDiffusion(Component):
 
     >>> mg = RasterModelGrid((5, 5), xy_spacing=2.0)
 
-    >>> c = mg.add_zeros("sediment_property__concentration", at="node")
-    >>> c[12] = 1.0
+    >>> mg.at_node["sediment_property__concentration"] = [
+    ...     [0.0, 0.0, 0.0, 0.0, 0.0],
+    ...     [0.0, 0.0, 0.0, 0.0, 0.0],
+    ...     [0.0, 0.0, 1.0, 0.0, 0.0],
+    ... ]
     >>> h = mg.add_full("soil__depth", 2.0, at="node")
     >>> z_br = mg.add_field(
     ...     "bedrock__elevation",
@@ -104,8 +110,9 @@ class ConcentrationTrackerForDiffusion(Component):
 
     >>> ddd = DepthDependentDiffuser(mg)
     >>> ct = ConcentrationTrackerForDiffusion(mg)
+    >>> ct.start_tracking()
     >>> ddd.run_one_step(1.0)
-    >>> ct.run_one_step(1.0)
+    >>> ct.stop_tracking(1.0)
 
     >>> mg.at_node["topographic__elevation"][mg.core_nodes].reshape((3, 3))
     array([[6. ,        7.13533528, 6. ],
@@ -118,8 +125,9 @@ class ConcentrationTrackerForDiffusion(Component):
 
     And running one more step.
 
+    >>> ct.start_tracking()
     >>> ddd.run_one_step(1.0)
-    >>> ct.run_one_step(1.0)
+    >>> ct.stop_tracking(1.0)
 
     >>> mg.at_node["topographic__elevation"][mg.core_nodes].reshape((3, 3))
     array([[5.52060315, 6.62473963, 5.52060315],
@@ -130,19 +138,22 @@ class ConcentrationTrackerForDiffusion(Component):
            [0.44750673, 1.        , 0.44750673],
            [0.09648071, 0.44750673, 0.09648071]])
 
-    Finally, the same 2D hillslope now using the DepthDependentTaylorDiffuser.
-    Note that the timestep must be smaller than 1 to maintain stability in the
-    diffusion calculation. Typically, one could use the dynamic timestepping
-    option. However, here it will provide incorrect soil flux values to the
-    ConcentrationTrackerForDiffusion, which cannot do sub-timestep calculations.
-    Use the if_unstable="warn" flag when instantiating the Taylor diffuser and
-    pick a timestep that is stable.
+    Finally, the same 2D hillslope now using the
+    :class:`~.DepthDependentTaylorDiffuser`. Note that the timestep must be smaller
+    than 1 to maintain stability in the diffusion calculation. Typically, one could
+    use the dynamic timestepping option. However, here it will provide incorrect
+    soil flux values to the :class:`~.ConcentrationTrackerForDiffusion`, which
+    cannot do sub-timestep calculations. Use the ``if_unstable="warn"`` flag when
+    instantiating the Taylor diffuser and pick a timestep that is stable.
 
     >>> from landlab.components import DepthDependentTaylorDiffuser
     >>> mg = RasterModelGrid((5, 5), xy_spacing=2.0)
 
-    >>> c = mg.add_zeros("sediment_property__concentration", at="node")
-    >>> c[12] = 1.0
+    >>> mg.at_node["sediment_property__concentration"] = [
+    ...     [0.0, 0.0, 0.0, 0.0, 0.0],
+    ...     [0.0, 0.0, 0.0, 0.0, 0.0],
+    ...     [0.0, 0.0, 1.0, 0.0, 0.0],
+    ... ]
     >>> h = mg.add_full("soil__depth", 2.0, at="node")
     >>> z_br = mg.add_field(
     ...     "bedrock__elevation",
@@ -154,8 +165,9 @@ class ConcentrationTrackerForDiffusion(Component):
 
     >>> ddtd = DepthDependentTaylorDiffuser(mg, if_unstable="warn")
     >>> ct = ConcentrationTrackerForDiffusion(mg)
+    >>> ct.start_tracking()
     >>> ddtd.run_one_step(0.4)
-    >>> ct.run_one_step(0.4)
+    >>> ct.stop_tracking(0.4)
 
     >>> mg.at_node["topographic__elevation"][mg.core_nodes].reshape((3, 3))
     array([[6.        , 7.30826823, 6.        ],
@@ -165,22 +177,13 @@ class ConcentrationTrackerForDiffusion(Component):
     array([[0.        , 0.26436925, 0.        ],
            [0.26436925, 1.        , 0.26436925],
            [0.        , 0.26436925, 0.        ]])
-
-    References
-    ----------
-    **Required Software Citation(s) Specific to this Component**
-
-    CITATION
-
     """
 
     _name = "ConcentrationTrackerForDiffusion"
 
     _unit_agnostic = True
 
-    _cite_as = """
-    CITATION
-    """
+    _cite_as = ""
 
     _info = {
         "soil__depth": {
@@ -254,7 +257,7 @@ def __init__(
         concentration_from_weathering: float or array, optional
             Concentration generated during the weathering process, -/m^3.
             Defaults to ``None``, which causes all weathered bedrock to retain its
-            original parent material concentration (``concentration_in_bedrock``)
+            original parent material concentration (`concentration_in_bedrock`)
             as it weathers to soil. Use this parameter to differentiate between
             the concentration in weathered material compared to its parent bedrock.
         """
@@ -328,7 +331,14 @@ def conc_w(self, new_val):
             raise ValueError("Concentration cannot be negative")
         self._conc_w = new_val
 
-    def calc_concentration(self, dt):
+    def _copy_old_soil_depth(self):
+        """Store a copy of soil depth. This is used as the old soil depth when
+        calculating changes in concentration.
+        """
+
+        self._soil__depth_old = self._soil__depth.copy()
+
+    def _calc_concentration(self, dt):
         """Calculate change in concentration for a time period 'dt'.
 
         Parameters
@@ -341,7 +351,6 @@ def calc_concentration(self, dt):
         conc_old = self._concentration.copy()
 
         # Map concentration from nodes to links (following soil flux direction)
-        # Does this overwrite fixed-value/gradient links?
         self._conc_at_links = map_value_at_max_node_to_link(
             self._grid, "topographic__elevation", "sediment_property__concentration"
         )
@@ -377,16 +386,28 @@ def calc_concentration(self, dt):
 
         self._concentration[~is_soil] = 0.0
 
-        # Update old soil depth to new value
-        self._soil__depth_old = self._soil__depth.copy()
+    def start_tracking(self):
+        """Stores values necessary for calculating changes in concentration.
+        This method must be called prior to running the sediment flux component
+        that changes physical properties of bedrock, soil, and/or topography.
+        """
 
-    def run_one_step(self, dt):
-        """Run for a time step.
+        self._copy_old_soil_depth()
+
+    def stop_tracking(self, dt):
+        """Calculate changes in concentration that have occurred in the timestep
+        since tracking was started. This method must be called after running the
+        sediment flux component that changes physical properties of bedrock,
+        soil, and/or topography.
 
         Parameters
         ----------
         dt: float (time)
             The imposed timestep.
         """
 
-        self.calc_concentration(dt)
+        self._calc_concentration(dt)
+
+    def run_one_step(self):
+        """run_one_step is not implemented for this component."""
+        raise NotImplementedError("run_one_step")

tests/components/concentration_tracker/test_concentration_tracker_for_diffusion.py

@@ -200,6 +200,13 @@ def setup_method(self):
         self.mg.add_zeros("soil_production__rate", at="node")
         self.mg.add_zeros("soil__flux", at="link")
 
+    def test_not_implemented(self):
+        """Test that private run_one_step is not implemented"""
+
+        ct = ConcentrationTrackerForDiffusion(self.mg)
+        with pytest.raises(NotImplementedError):
+            ct.run_one_step()
+
     def test_concentration_from_soil_flux(self):
         """
         ConcentrationTrackerForDiffusion should correctly calculate concentration
@@ -220,7 +227,8 @@ def test_concentration_from_soil_flux(self):
         # Flux of -1 in the middle row should shift all C_sed values left by one node.
 
         ct = ConcentrationTrackerForDiffusion(self.mg)
-        ct.run_one_step(1)
+        ct.start_tracking()
+        ct.stop_tracking(1)
 
         expected = np.asarray(
             [
@@ -256,14 +264,13 @@ def test_concentration_from_weathering_without_conc_w(self):
 
         ct = ConcentrationTrackerForDiffusion(self.mg)
 
+        ct.start_tracking()
         # Soil volume is 1 at each node. Soil production rate of 1 doubles volume.
         # This is normally done by the DepthDependentDiffuser. Here, it is forced.
         self.mg.at_node["soil__depth"] += 1.0
-
         # Node 7: C_sed remains 1 because parent bedrock had conc_br of 1.
         # Node 8: C_sed is halved from 1 to 0.5 because parent bedrock had conc_br = 0.
-
-        ct.run_one_step(1)
+        ct.stop_tracking(1)
 
         # Node 7 should have the same concentration as before.
         # Node 8 should have half its previous concentration.
@@ -301,17 +308,16 @@ def test_concentration_from_weathering_with_conc_w(self):
             self.mg, concentration_from_weathering=0.0
         )
 
+        ct.start_tracking()
         # Soil volume is 1 at each node. Soil production rate of 1 doubles volume.
         # This is normally done by the DepthDependentDiffuser. Here, it is forced.
         self.mg.at_node["soil__depth"] += 1.0
-
         # conc_w overrides conc_br values. In this case, no concentration is produced by
         # the weathering process, even at Node 7 where conc_br = 1.
 
         # Node 7: C_sed is halved from 1 to 0.5 despite parent bedrock with conc_br = 1.
         # Node 8: C_sed is halved from 1 to 0.5 because conc_w = 0.
-
-        ct.run_one_step(1)
+        ct.stop_tracking(1)
 
         # Node 7 should have half its previous concentration.
         # Node 8 should have half its previous concentration.
@@ -373,7 +379,8 @@ def test_mass_balance_with_all_boundaries_closed(self):
         )
 
         ct = ConcentrationTrackerForDiffusion(self.mg)
-        ct.run_one_step(1)
+        ct.start_tracking()
+        ct.stop_tracking(1)
 
         total_mass_after = np.sum(
             self.mg.at_node["sediment_property__concentration"]
@@ -397,7 +404,8 @@ def test_mass_balance_with_one_boundary_open(self):
         )
 
         ct = ConcentrationTrackerForDiffusion(self.mg)
-        ct.run_one_step(1)
+        ct.start_tracking()
+        ct.stop_tracking(1)
 
         total_mass_leaving = 1
         total_mass_after = np.sum(
@@ -407,3 +415,33 @@ def test_mass_balance_with_one_boundary_open(self):
         )
 
         assert_allclose(total_mass_before, total_mass_after + total_mass_leaving)
+
+
+# %%
+class TestFieldCopy:
+    """Test that copied field is a copy, but not a reference."""
+
+    def setup_method(self):
+        self.mg = RasterModelGrid((3, 5))
+        self.mg.add_zeros("soil__flux", at="link")
+        self.mg.add_zeros("soil__depth", at="node")
+        self.mg.add_zeros("soil_production__rate", at="node")
+        self.mg.add_zeros("topographic__elevation", at="node")
+
+    def test_copy_is_equal(self):
+        """Test that copied values are equal to copied field."""
+
+        ct = ConcentrationTrackerForDiffusion(self.mg)
+        ct._copy_old_soil_depth()
+
+        assert np.allclose(ct._soil__depth_old, self.mg.at_node["soil__depth"])
+
+    def test_copy_is_not_reference(self):
+        """Test that copy not a reference."""
+
+        ct = ConcentrationTrackerForDiffusion(self.mg)
+        ct._copy_old_soil_depth()
+
+        self.mg.at_node["soil__depth"] += 1
+
+        assert not np.allclose(ct._soil__depth_old, self.mg.at_node["soil__depth"])