Source code for landlab.components.gflex.flexure

"""This is a Landlab wrapper for A Wickert's gFlex flexure model (Wickert et
al., submitted to Geoscientific Model Development). The most up-to-date version
of his code can be found at github.com/awickert/gFlex.

This Landlab wrapper will use a snapshot of that code, which YOU need to
install on your own machine.
A stable snapshot of gFlex is hosted on PyPI, which is the recommended version
to install.
If you have pip (the Python package install tool), simply run
'pip install gFlex' from a command prompt.
Alternatively, you can download and unpack the code (from github, or with PyPI,
pypi.python.org/pypi/gFlex/), then run 'python setup.py install'.

Created on Thu Feb 19 18:47:11 2015

@author: daniel.hobley (SiccarPoint @Github)

...following AW's run_in_script_2D.py.
"""

import numpy as np
import scipy.constants

from landlab import Component
from landlab import FieldError
from landlab import RasterModelGrid

try:
    import gflex
except ImportError:
    NO_GFLEX = True
else:
    NO_GFLEX = False


[docs] class gFlex(Component): """This is a Landlab wrapper for A Wickert's gFlex flexure model (Wickert et al., 2016, Geoscientific Model Development). The most up-to-date version of his code can be found at github.com/awickert/gFlex. This Landlab wrapper will use a snapshot of that code, which YOU need to install on your own machine. A stable snapshot of gFlex is hosted on PyPI, which is the recommended version to install. If you have pip (the Python package install tool), simply run 'pip install gFlex' from a command prompt. Alternatively, you can download and unpack the code (from github, or with PyPI, pypi.python.org/pypi/gFlex/), then run 'python setup.py install'. Note that gFlex maintains its own internal version if the grid, but this should not affect performance. This component will modify the topographic__elevation field only if one already exists. Note that the gFlex component **demands lengths in meters**, including the grid dimensions. The component also recognises the gFlex specific parameters 'Method', 'PlateSolutionType', 'Solver', and 'Quiet'. See the gFlex software documentation for more details. Examples -------- NB: these tests are not actually run as our automated testing becomes confused if gFlex is not installed on the testing machine! >>> from landlab import RasterModelGrid >>> from landlab.components import gFlex >>> mg = RasterModelGrid((10, 10), xy_spacing=25000.0) >>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float) >>> stress = mg.add_zeros("surface_load__stress", at="node", dtype=float) >>> stress.view().reshape(mg.shape)[3:7, 3:7] += 1.0e6 >>> gf = gFlex( ... mg, BC_E="0Moment0Shear", BC_N="Periodic", BC_S="Periodic" ... ) # doctest: +SKIP >>> gf.run_one_step() # doctest: +SKIP N-S profile across flexed plate: >>> z.reshape(mg.shape)[:, 5] # doctest: +SKIP array([-4.54872677, -4.6484927 , -4.82638669, -5.03001546, -5.15351385, -5.15351385, -5.03001546, -4.82638669, -4.6484927 , -4.54872677]) W-E profile, noting the free BC to the east side: >>> z.reshape(mg.shape)[5, :] # doctest: +SKIP array([-0.43536739, -1.19197738, -2.164915 , -3.2388464 , -4.2607558 , -5.15351385, -5.89373366, -6.50676947, -7.07880156, -7.63302576]) References ---------- **Required Software Citation(s) Specific to this Component** Wickert, A. (2016). Open-source modular solutions for flexural isostasy: gFlex v1.0. Geoscientific Model Development 9(3), 997-1017. https://dx.doi.org/10.5194/gmd-9-997-2016 **Additional References** None Listed """ _name = "gFlex" _unit_agnostic = True _cite_as = """ @article{wickert2016open, author = {Wickert, A. D.}, title = {{Open-source modular solutions for flexural isostasy: gFlex v1.0}}, issn = {1991-959X}, doi = {10.5194/gmd-9-997-2016}, pages = {997--1017}, number = {3}, volume = {9}, journal = {Geoscientific Model Development}, year = {2016} } """ _info = { "lithosphere_surface__elevation_increment": { "dtype": float, "intent": "out", "optional": False, "units": "m", "mapping": "node", "doc": ( "The change in elevation of the top of the lithosphere (the " "land surface) in one timestep" ), }, "surface_load__stress": { "dtype": float, "intent": "in", "optional": False, "units": "Pa", "mapping": "node", "doc": "Magnitude of stress exerted by surface load", }, "topographic__elevation": { "dtype": float, "intent": "out", "optional": False, "units": "m", "mapping": "node", "doc": "Land surface topographic elevation", }, }
[docs] def __init__( self, grid, Youngs_modulus=6.5e11, Poissons_ratio=0.25, rho_mantle=3300.0, rho_fill=0.0, elastic_thickness=35000.0, Method="FD", Solver="direct", PlateSolutionType="vWC1994", quiet=True, BC_W="0Displacement0Slope", BC_E="0Displacement0Slope", BC_N="0Displacement0Slope", BC_S="0Displacement0Slope", g=scipy.constants.g, ): """Constructor for Wickert's gFlex in Landlab. Parameters ---------- Youngs_modulus : float Young's modulus for the lithosphere. Poissons_ratio : float Poisson's ratio for the lithosphere. rho_mantle : float (kg*m**-3) The density of the mantle. rho_fill : float (kg*m**-3) The density of the infilling material (air, water...) elastic_thickness : float (m) The elastic thickness of the lithosphere. BC_W, BC_E, BC_N, BC_S : {'0Displacement0Slope', '0Moment0Shear', 'Periodic'} The boundary condition status of each grid edge, following gFlex's definitions. Periodic boundaries must be paired (obviously). g : float (m*s**-2) The acceleration due to gravity. """ super().__init__(grid) assert isinstance(grid, RasterModelGrid) if NO_GFLEX: raise ImportError( "gFlex not installed! For installation instructions see " + "gFlex on GitHub: https://github.com/awickert/gFlex" ) BC_options = ( "0Displacement0Slope", "0Moment0Shear", "0Slope0Shear", "Periodic", ) # instantiate the module: self._flex = gflex.F2D() flex = self._flex # set up the grid variables: flex.dx = grid.dx flex.dy = grid.dy # we assume these properties are fixed in this relatively # straightforward implementation, but they can still be set if you # want: flex.Method = Method flex.PlateSolutionType = PlateSolutionType flex.Solver = Solver flex.Quiet = quiet flex.E = float(Youngs_modulus) flex.nu = float(Poissons_ratio) flex.rho_m = float(rho_mantle) flex.rho_fill = float(rho_fill) flex.g = float(g) flex.BC_W = BC_W flex.BC_E = BC_E flex.BC_S = BC_S flex.BC_N = BC_N for i in (flex.BC_E, flex.BC_W, flex.BC_N, flex.BC_S): assert i in BC_options if BC_W == "Periodic": assert BC_E == "Periodic" if BC_E == "Periodic": assert BC_W == "Periodic" if BC_N == "Periodic": assert BC_S == "Periodic" if BC_S == "Periodic": assert BC_N == "Periodic" Te_in = elastic_thickness try: flex.Te = float(Te_in) except ValueError: try: flex.Te = grid.at_node[Te_in].view().reshape(grid.shape) except TypeError: flex.Te = Te_in.view().reshape(grid.shape) self._input_var_names.add(Te_in) self._output_var_names.add(Te_in) # set up the link between surface load stresses in the gFlex component # and the LL grid field: flex.qs = grid.at_node["surface_load__stress"].view().reshape(grid.shape) # create a holder for the "pre-flexure" state of the grid, to allow # updating of elevs: self._pre_flex = np.zeros(grid.number_of_nodes, dtype=float) # create the primary output field: self._grid.add_zeros( "lithosphere_surface__elevation_increment", at="node", dtype=float, clobber=True, )
[docs] def flex_lithosphere(self): """Executes (& finalizes, from the perspective of gFlex) the core method of gFlex. Note that flexure of the lithosphere proceeds to steady state in a single timestep. """ self._flex.qs = ( self._grid.at_node["surface_load__stress"].view().reshape(self._grid.shape) ) self._flex.initialize() self._flex.run() self._flex.finalize() self._grid.at_node["lithosphere_surface__elevation_increment"][ : ] = self._flex.w.view().ravel() try: self._grid.at_node["topographic__elevation"] # a topo exists... except FieldError: pass else: topo_diff = ( self._grid.at_node["lithosphere_surface__elevation_increment"] - self._pre_flex ) self._grid.at_node["topographic__elevation"] += topo_diff self._pre_flex += topo_diff
[docs] def run_one_step(self): """Flex the lithosphere to find its steady state form. The standardized run method for this component. Parameters ---------- None """ self._flex_lithosphere()