landlab.components.diffusion.diffusion

Component that models 2D diffusion using an explicit finite-volume method.

Created July 2013 GT Last updated March 2016 DEJH with LL v1.0 component style

class LinearDiffuser

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)

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.

static __new__(cls, *args, **kwds)

cite_as = ''

property coords

Return the coordinates of nodes on grid attached to the component.

property current_time

Current time.

Some components may keep track of the current time. In this case, the current_time attribute is incremented. Otherwise it is set to None.

Return type:

current_time

definitions = (('hillslope_sediment__unit_volume_flux', 'Volume flux per unit width along links'), ('topographic__elevation', 'Land surface topographic elevation'), ('topographic__gradient', 'Gradient of the ground surface'))

property fixed_grad_anchors

Fixed gradient anchors.

property fixed_grad_nodes

Fixed gradient nodes.

property fixed_grad_offsets

Fixed gradient offsets.

classmethod from_path(grid, path)

Create a component from an input file.

Parameters:

• grid (ModelGrid) – A landlab grid.

• path (str or file_like) – Path to a parameter file, contents of a parameter file, or a file-like object.

Returns:

A newly-created component.

Return type:

Component

property grid

Return the grid attached to the component.

initialize_optional_output_fields()

Create fields for a component based on its optional field outputs, if declared in _optional_var_names.

This method will create new fields (without overwrite) for any fields output by the component as optional. New fields are initialized to zero. New fields are created as arrays of floats, unless the component also contains the specifying property _var_type.

initialize_output_fields(values_per_element=None)

Create fields for a component based on its input and output var names.

This method will create new fields (without overwrite) for any fields output by, but not supplied to, the component. New fields are initialized to zero. Ignores optional fields. New fields are created as arrays of floats, unless the component specifies the variable type.

Parameters:

values_per_element (int (optional)) – On occasion, it is necessary to create a field that is of size (n_grid_elements, values_per_element) instead of the default size (n_grid_elements,). Use this keyword argument to acomplish this task.

input_var_names = ('topographic__elevation',)

name = 'LinearDiffuser'

optional_var_names = ()

output_var_names = ('hillslope_sediment__unit_volume_flux', 'topographic__elevation', 'topographic__gradient')

run_one_step(dt)

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 shape

Return the grid shape attached to the component, if defined.

property time_step

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

unit_agnostic = True

units = (('hillslope_sediment__unit_volume_flux', 'm**2/s'), ('topographic__elevation', 'm'), ('topographic__gradient', '-'))

updated_boundary_conditions()

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

classmethod var_definition(name)

Get a description of a particular field.

Parameters:

name (str) – A field name.

Returns:

A description of each field.

Return type:

tuple of (name, *description*)

classmethod var_help(name)

Print a help message for a particular field.

Parameters:

name (str) – A field name.

classmethod var_loc(name)

Location where a particular variable is defined.

Parameters:

name (str) – A field name.

Returns:

The location (‘node’, ‘link’, etc.) where a variable is defined.

Return type:

str

var_mapping = (('hillslope_sediment__unit_volume_flux', 'link'), ('topographic__elevation', 'node'), ('topographic__gradient', 'link'))

classmethod var_type(name)

Returns the dtype of a field (float, int, bool, str…).

Parameters:

name (str) – A field name.

Returns:

The dtype of the field.

Return type:

dtype

classmethod var_units(name)

Get the units of a particular field.

Parameters:

name (str) – A field name.

Returns:

Units for the given field.

Return type:

str