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

class LinearDiffuser(self, grid, linear_diffusivity=None, method='simple', deposit=True, **kwds)[source]

Bases: landlab.core.model_component.Component

This component implements 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. Options other than ‘simple’ will make the component run slower, but give second order solutions that incorporate information from more distal nodes (i.e., the diagonals). Note that the option ‘resolve_on_patches’ can result in somewhat counterintuitive behaviour - 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.


>>> from landlab import RasterModelGrid
>>> import numpy as np
>>> mg = RasterModelGrid((9, 9))
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> z.reshape((9, 9))[4, 4] = 1.
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, True)
>>> ld = LinearDiffuser(mg, linear_diffusivity=1.)
>>> for i in range(1):
...     ld.run_one_step(1.)
>>> np.isclose(z[mg.core_nodes].sum(), 1.)
>>> mg2 = RasterModelGrid((5, 30))
>>> z2 = mg2.add_zeros('node', 'topographic__elevation')
>>> z2.reshape((5, 30))[2, 8] = 1.
>>> z2.reshape((5, 30))[2, 22] = 1.
>>> mg2.set_closed_boundaries_at_grid_edges(True, True, True, True)
>>> kd = mg2.node_x/mg2.node_x.mean()
>>> ld2 = LinearDiffuser(mg2, linear_diffusivity=kd)
>>> for i in range(10):
...     ld2.run_one_step(0.1)
>>> z2[mg2.core_nodes].sum() == 2.
>>> z2.reshape((5, 30))[2, 8] > z2.reshape((5, 30))[2, 22]

An example using links:

>>> mg1 = RasterModelGrid((10, 10), xy_spacing=100.)
>>> mg2 = RasterModelGrid((10, 10), xy_spacing=100.)
>>> z1 = mg1.add_zeros('node', 'topographic__elevation')
>>> z2 = mg2.add_zeros('node', 'topographic__elevation')
>>> dt = 1.
>>> nt = 10
>>> kappa_links = mg2.add_ones('link', 'surface_water__discharge')
>>> kappa_links *= 10000.
>>> dfn1 = LinearDiffuser(mg1, linear_diffusivity=10000.)
>>> dfn2 = LinearDiffuser(mg2, linear_diffusivity='surface_water__discharge')
>>> for i in range(nt):
...     z1[mg1.core_nodes] += 1.
...     z2[mg2.core_nodes] += 1.
...     dfn1.run_one_step(dt)
...     dfn2.run_one_step(dt)
>>> np.allclose(z1, z2)
>>> z2.fill(0.)
>>> dfn2 = LinearDiffuser(mg2, linear_diffusivity='surface_water__discharge',
...                       method='resolve_on_patches')
>>> for i in range(nt):
...     z2[mg2.core_nodes] += 1.
...     dfn2.run_one_step(dt)
>>> np.all(z2[mg2.core_nodes] < z1[mg2.core_nodes])
  • 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', 'on_diagonals'}) – 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. ‘on_diagonals’ permits Raster diagonals to carry the diffusional fluxes. These latter techniques are more computationally expensive, but can suppress cardinal direction artifacts if diffusion is performed on a raster. ‘on_diagonals’ pretends that the “faces” of a cell with 8 links are represented by a stretched regular octagon set within the true cell.

  • 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.

diffuse(dt, **kwds)[source]

See run_one_step.

run_one_step(dt, **kwds)[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.


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

property time_step

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


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


>>> from landlab import RasterModelGrid
>>> import numpy as np
>>> mg = RasterModelGrid((4, 5))
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> z[mg.core_nodes] = 1.
>>> ld = LinearDiffuser(mg, linear_diffusivity=1.)
>>> ld.fixed_grad_nodes.size == 0
>>> ld.fixed_grad_anchors.size == 0
>>> ld.fixed_grad_offsets.size == 0
>>> mg.at_link['topographic__slope'] = mg.calc_grad_at_link(
...     'topographic__elevation')
>>> mg.set_fixed_link_boundaries_at_grid_edges(True, True, True, True)
>>> 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)