In this example we will:
A bit of magic so that we can plot within this notebook.
In [ ]:
%matplotlib inline
import numpy as np
We are going to build a uniform rectilinear grid with a node spacing of 10 km in the y-direction and 20 km in the x-direction on which we will solve the flexure equation.
First we need to import RasterModelGrid
. We also import the Landlab plotting function imshow_grid
to view the grids.
In [ ]:
from landlab import RasterModelGrid
from landlab.plot.imshow import imshow_grid
Create a rectilinear grid with a spacing of 10 km between rows and 20 km between columns. The numbers of rows and columms are provided as a tuple
of (n_rows, n_cols)
, in the same manner as similar numpy functions. The spacing is also a tuple
, (dy, dx)
.
In [ ]:
grid = RasterModelGrid((200, 400), xy_spacing=(10e3, 20e3))
In [ ]:
grid.dy, grid.dx
Now we create the flexure component and tell it to use our newly-created grid. First, though, we'll examine the Flexure
component a bit.
In [ ]:
from landlab.components.flexure import Flexure
The Flexure component, as with most landlab components, will require our grid to have some data that it will use. We can get the names of these data fields with the intput_var_names
attribute of the component class.
In [ ]:
Flexure.input_var_names
We see that flexure uses just one data field: the change in lithospheric loading. landlab component classes can provide additional information about each of these fields. For instance, to the the units for a field, use the var_units
method.
In [ ]:
Flexure.var_units('lithosphere__overlying_pressure_increment')
To print a more detailed description of a field, use var_help
.
In [ ]:
Flexure.var_help('lithosphere__overlying_pressure_increment')
What about the data that Flexure
provides? Use the output_var_names
attribute.
In [ ]:
Flexure.output_var_names
In [ ]:
Flexure.var_help('lithosphere_surface__elevation_increment')
Now that we understand the component a little more, create it using our grid.
In [ ]:
grid.add_zeros("lithosphere__overlying_pressure_increment", at="node")
flex = Flexure(grid, method='flexure', n_procs=4)
In [ ]:
load = np.random.normal(0, 100 * 2650. * 9.81, grid.number_of_nodes)
grid.at_node['lithosphere__overlying_pressure_increment'] = load
In [ ]:
imshow_grid(grid,
'lithosphere__overlying_pressure_increment',
symmetric_cbar=True,
cmap='nipy_spectral')
In [ ]:
flex.update()
As we saw above, the flexure component creates an output field (lithosphere_surface__elevation_increment
) that contains surface deflections for the applied loads.
We now plot these deflections with the imshow_grid
method, which is available to all landlab components.
In [ ]:
imshow_grid(grid,
'lithosphere_surface__elevation_increment',
symmetric_cbar=True,
cmap='nipy_spectral')
Maintain the same loading distribution but double the effective elastic thickness.
In [ ]:
flex.eet *= 2.
flex.update()
imshow_grid(grid,
'lithosphere_surface__elevation_increment',
symmetric_cbar=True,
cmap='nipy_spectral')
Now let's add a vertical rectangular load to the middle of the grid. We plot the load grid first to make sure we did this correctly.
In [ ]:
load[np.where(np.logical_and(grid.node_x>3000000, grid.node_x<5000000))]= \
load[np.where(np.logical_and(grid.node_x>3000000, grid.node_x<5000000))]+1e7
imshow_grid(grid,
'lithosphere__overlying_pressure_increment',
symmetric_cbar=True,
cmap='nipy_spectral')
In [ ]:
flex.update()
imshow_grid(grid,
'lithosphere_surface__elevation_increment',
symmetric_cbar=True,
cmap='nipy_spectral')
In [ ]: