LinearDiffuser: Model soil creep using “linear diffusion” transport law (no depth dependence)#

class LinearDiffuser(*args, **kwds)[source]#

Bases: Component

Linear diffusion of a Landlab field.

Component assumes grid does not deform. If the boundary conditions on the grid change after component instantiation, be sure to also call updated_boundary_conditions to ensure these are reflected in the component (especially if fixed_links are present).

The method keyword allows control of the way the solver works. Note that the option ‘resolve_on_patches’ can result in somewhat counterintuitive behavior - this option tells the component to treat the diffusivity as a field with directionality to it (i.e., that the diffusivites are defined on links). Thus if all links have the same diffusivity value, with this flag active “effective” diffusivities at the nodes will be higher than this value (by a factor of root 2) as the diffusivity at each patch will be the mean vector sum of that at the bounding links.

The primary method of this class is run_one_step.

Examples

>>> from landlab import RasterModelGrid
>>> import numpy as np

>>> grid = RasterModelGrid((9, 9))
>>> z = grid.add_zeros("topographic__elevation", at="node")
>>> z.reshape((9, 9))[4, 4] = 1.0
>>> grid.set_closed_boundaries_at_grid_edges(True, True, True, True)
>>> ld = LinearDiffuser(grid, linear_diffusivity=1.0)
>>> ld.run_one_step(1.0)
>>> np.isclose(z[grid.core_nodes].sum(), 1.0)
True

>>> grid = RasterModelGrid((5, 30))
>>> z2 = grid.add_zeros("topographic__elevation", at="node")
>>> z2.reshape((5, 30))[2, 8] = 1.0
>>> z2.reshape((5, 30))[2, 22] = 1.0
>>> grid.set_closed_boundaries_at_grid_edges(True, True, True, True)
>>> kd = grid.node_x / grid.node_x.mean()
>>> ld2 = LinearDiffuser(grid, linear_diffusivity=kd)
>>> for _ in range(10):
...     ld2.run_one_step(0.1)
...
>>> z2[grid.core_nodes].sum() == 2.0
True
>>> z2.reshape((5, 30))[2, 8] > z2.reshape((5, 30))[2, 22]
True

An example using links:

>>> grid1 = RasterModelGrid((10, 10), xy_spacing=100.0)
>>> grid2 = RasterModelGrid((10, 10), xy_spacing=100.0)
>>> z1 = grid1.add_zeros("topographic__elevation", at="node")
>>> z2 = grid2.add_zeros("topographic__elevation", at="node")
>>> dt = 1.0
>>> nt = 10
>>> grid2.at_link["surface_water__discharge"] = np.full(
...     grid2.number_of_links, 10000.0
... )
>>> dfn1 = LinearDiffuser(grid1, linear_diffusivity=10000.0)
>>> dfn2 = LinearDiffuser(grid2, linear_diffusivity="surface_water__discharge")
>>> for i in range(nt):
...     z1[grid1.core_nodes] += 1.0
...     z2[grid2.core_nodes] += 1.0
...     dfn1.run_one_step(dt)
...     dfn2.run_one_step(dt)
...
>>> np.allclose(z1, z2)
True
>>> z2.fill(0.0)
>>> dfn2 = LinearDiffuser(
...     grid2,
...     linear_diffusivity="surface_water__discharge",
...     method="resolve_on_patches",
... )
>>> for i in range(nt):
...     z2[grid2.core_nodes] += 1.0
...     dfn2.run_one_step(dt)
...
>>> np.all(z2[grid2.core_nodes] < z1[grid2.core_nodes])
True

References

Required Software Citation(s) Specific to this Component

None Listed

Additional References

Culling, W. (1963). Soil Creep and the Development of Hillside Slopes. The Journal of Geology 71(2), 127-161. https://dx.doi.org/10.1086/626891

Parameters:

  • grid (ModelGrid) – A grid.

  • linear_diffusivity (float, array, or field name (m**2/time)) – The diffusivity. If an array or field name, these must be the diffusivities on either nodes or links - the component will distinguish which based on array length. Values on nodes will be mapped to links using an upwind scheme in the simple case.

  • method ({'simple', 'resolve_on_patches'}) – The method used to represent the fluxes. ‘simple’ solves a finite difference method with a simple staggered grid scheme onto the links. ‘resolve_on_patches’ solves the scheme by mapping both slopes and diffusivities onto the patches and solving there before resolving values back to the nodes (and at the moment requires a raster grid). Note that this scheme enforces directionality in the diffusion field; it’s no longer just a scalar field. Thus diffusivities must be defined on links when this option is chosen.

  • deposit ({True, False}) – Whether diffusive material can be deposited. True means that diffusive material will be deposited if the divergence of sediment flux is negative. False means that even when the divergence of sediment flux is negative, no material is deposited. (No deposition ever.) The False case is a bit of a band-aid to account for cases when fluvial incision likely removes any material that would be deposited. If one couples fluvial detachment-limited incision with linear diffusion, the channels will not reach the predicted analytical solution unless deposit is set to False.

__init__(grid, linear_diffusivity=0.01, method='simple', deposit=True)[source]#
Parameters:

  • grid (ModelGrid) – A grid.

  • linear_diffusivity (float, array, or field name (m**2/time)) – The diffusivity. If an array or field name, these must be the diffusivities on either nodes or links - the component will distinguish which based on array length. Values on nodes will be mapped to links using an upwind scheme in the simple case.

  • method ({'simple', 'resolve_on_patches'}) – The method used to represent the fluxes. ‘simple’ solves a finite difference method with a simple staggered grid scheme onto the links. ‘resolve_on_patches’ solves the scheme by mapping both slopes and diffusivities onto the patches and solving there before resolving values back to the nodes (and at the moment requires a raster grid). Note that this scheme enforces directionality in the diffusion field; it’s no longer just a scalar field. Thus diffusivities must be defined on links when this option is chosen.

  • deposit ({True, False}) – Whether diffusive material can be deposited. True means that diffusive material will be deposited if the divergence of sediment flux is negative. False means that even when the divergence of sediment flux is negative, no material is deposited. (No deposition ever.) The False case is a bit of a band-aid to account for cases when fluvial incision likely removes any material that would be deposited. If one couples fluvial detachment-limited incision with linear diffusion, the channels will not reach the predicted analytical solution unless deposit is set to False.

property fixed_grad_anchors#

Fixed gradient anchors.

property fixed_grad_nodes#

Fixed gradient nodes.

property fixed_grad_offsets#

Fixed gradient offsets.

run_one_step(dt)[source]#

Run the diffuser for one timestep, dt.

If the imposed timestep dt is longer than the Courant-Friedrichs-Lewy condition for the diffusion, this timestep will be internally divided as the component runs, as needed.

Parameters:

dt (float (time)) – The imposed timestep.

property time_step#

Returns internal time-step size (as a property).

updated_boundary_conditions()[source]#

Call if grid BCs are updated after component instantiation.

Sets fixed_grad_nodes, fixed_grad_anchors, & fixed_grad_offsets, such that:

value[fixed_grad_nodes] = value[fixed_grad_anchors] + offset

Examples

>>> from landlab import RasterModelGrid
>>> import numpy as np
>>> mg = RasterModelGrid((4, 5))
>>> z = mg.add_zeros("topographic__elevation", at="node")
>>> z[mg.core_nodes] = 1.0
>>> ld = LinearDiffuser(mg, linear_diffusivity=1.0)
>>> ld.fixed_grad_nodes.size == 0
True
>>> ld.fixed_grad_anchors.size == 0
True
>>> ld.fixed_grad_offsets.size == 0
True
>>> mg.at_link["topographic__slope"] = mg.calc_grad_at_link(
...     "topographic__elevation"
... )
>>> mg.status_at_node[mg.perimeter_nodes] = mg.BC_NODE_IS_FIXED_GRADIENT
>>> ld.updated_boundary_conditions()
>>> ld.fixed_grad_nodes
array([ 1,  2,  3,  5,  9, 10, 14, 16, 17, 18])
>>> ld.fixed_grad_anchors
array([ 6,  7,  8,  6,  8, 11, 13, 11, 12, 13])
>>> ld.fixed_grad_offsets
array([-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.])
>>> np.allclose(
...     z[ld.fixed_grad_nodes], z[ld.fixed_grad_anchors] + ld.fixed_grad_offsets
... )
True