making_components/making_components.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"http://landlab.github.io\"><img style=\"float: left; height: 175px; width: 175px\" src=\"../landlab_logo_picture.jpg\"></a> <h3 style=\"margin: 117px 0 0 185px; font-weight: 300;\">a toolkit for modeling earth surface processes</h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to write a Landlab component"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "<small> For instructions on how to run an interactive iPython notebook, click here: <a href=\"https://github.com/landlab/tutorials/blob/release/README.md\">https://github.com/landlab/tutorials/blob/release/README.md</a></small><br>\n",
    "<small>For the unexpanded version to download and run, click here: <a href=\"http://nbviewer.jupyter.org/github/landlab/tutorials/blob/release/making_components/making_components_unexpanded.ipynb\">http://nbviewer.jupyter.org/github/landlab/tutorials/blob/release/making_components/making_components_unexpanded.ipynb</a></small><br>\n",
    "<small>For more Landlab tutorials, click here: <a href=\"https://github.com/landlab/landlab/wiki/Tutorials\">https://github.com/landlab/landlab/wiki/Tutorials</a></small>\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [NOTE: THIS TUTORIAL PERTAINS TO LANDLAB 1.x. SOME ASPECTS OF COMPONENT DESIGN, SUCH AS HEADER INFORMATION, WILL BE RE-DESIGNED IN LANDLAB 2.0]\n",
    "\n",
    "This ipython notebook walks you through the basic procedure for writing a Landlab component, using the example of a kinematic-wave flow model.\n",
    "\n",
    "## Overview\n",
    "A Landlab component is implemented as a Python `class`. Although every Landlab component is unique in some respects, to be a component, a class must have at least the following standard ingredients:\n",
    "\n",
    "(1) The class must inherit the base class `Component`.\n",
    "\n",
    "(2) The class must include a set of standard variables defined in the header (i.e., before the `__init__` method), which describe the data arrays that the component uses.\n",
    "\n",
    "(3) The class must have an `__init__` method defined, with a semi-standardized parameter list described below.\n",
    "\n",
    "(4) The class must provide a function that does performs the component's \"action\", typically named `run_one_step()` and this function's parameter list must follow the convention described below.\n",
    "\n",
    "\n",
    "## Class definition and header\n",
    "\n",
    "A Landlab component is a class that inherits from `Component`. The name of the class should be in CamelCase, and should make sense when used in the sentence: \"A *(component-name)* is a...\". The class definition should be followed by a docstring. The docstring should include a list of parameters for the `__init__` method and succintly describe them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from landlab import Component, FieldError\n",
    "import numpy as np\n",
    "\n",
    "class KinwaveOverlandFlowModel(Component):\n",
    "    \"\"\"\n",
    "    Calculate water flow over topography.\n",
    "    \n",
    "    Landlab component that implements a two-dimensional \n",
    "    kinematic wave model.\n",
    "    \n",
    "    Construction:\n",
    "    \n",
    "        KinwaveOverlandFlowModel(grid, precip_rate=1.0,\n",
    "                                 precip_duration=1.0,\n",
    "                                 infilt_rate=0.0,\n",
    "                                 roughness=0.01)\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    grid : ModelGrid\n",
    "        A Landlab grid object.\n",
    "    precip_rate : float, optional (defaults to 1 mm/hr)\n",
    "        Precipitation rate, mm/hr\n",
    "    precip_duration : float, optional (defaults to 1 hour)\n",
    "        Duration of precipitation, hours\n",
    "    infilt_rate : float, optional (defaults to 0)\n",
    "        Maximum rate of infiltration, mm/hr\n",
    "    roughnes : float, defaults to 0.01\n",
    "        Manning roughness coefficient, s/m^1/3\n",
    "\n",
    "    \"\"\"\n",
    "    def __init__(): # ignore this for now, we will add more stuff eventually. \n",
    "        pass "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Doc tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following docstring section 'Examples' should help the user understand what is the component's purpose and how it works. It is an example (or examples) of its use in a (more or less) simple case within the Landlab framework: a grid is created, the component is instantiated on this grid and run. Unlike in the example below, we strongly recommend commenting your example(s) to explain what is happening.\n",
    "\n",
    "This is also the section that will be run during the integration tests of your component (once you have submitted a pull request to have your component merged into the Landlab release branch). All lines starting with >>> are run and should produce the results you provided: here, the test will fail if `kw.vel_coeff` does not return `100.0.`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\nExamples\\n--------\\n>>> from landlab import RasterModelGrid\\n>>> rg = RasterModelGrid((4, 5), 10.0)\\n>>> kw = KinwaveOverlandFlowModel(rg)\\n>>> kw.vel_coef\\n100.0\\n>>> rg.at_node['surface_water__depth']\\narray([ 0.,  0.,  0.,  0.,  0.,\\n        0.,  0.,  0.,  0.,  0.,\\n        0.,  0.,  0.,  0.,  0.,\\n        0.,  0.,  0.,  0.,  0.])\\n\""
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "    \"\"\"\n",
    "    Examples\n",
    "    --------\n",
    "    >>> from landlab import RasterModelGrid\n",
    "    >>> rg = RasterModelGrid((4, 5), 10.0)\n",
    "    >>> kw = KinwaveOverlandFlowModel(rg)\n",
    "    >>> kw.vel_coef\n",
    "    100.0\n",
    "    >>> rg.at_node['surface_water__depth']\n",
    "    array([ 0.,  0.,  0.,  0.,  0.,\n",
    "            0.,  0.,  0.,  0.,  0.,\n",
    "            0.,  0.,  0.,  0.,  0.,\n",
    "            0.,  0.,  0.,  0.,  0.])\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Header information: `_name`\n",
    "Every component should have a name, as a string. Normally this will be the same as the class name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "_name = 'KinwaveOverlandFlowModel'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Header information: `_input_var_names`\n",
    "\n",
    "(NOTE: IN LANDLAB 2.0, THIS AND OTHER HEADER VARIABLES BELOW WILL BE REPLACED WITH A VARIABLE CALLED `_info` THAT CONTAINS A DICTIONARY OF FIELD NAMES. FOR EACH ENTRY, THERE IS A SUB-DICTIONARY CONTAINING DATA PERTAINING TO THAT FIELD, INCLUDING ITS INTENT (INPUT OR OUTPUT), ITS UNITS, ITS MAPPING, ETC.) \n",
    "\n",
    "\"Input variables\" here refers to those arrays (fields) that the component needs, and for which initial values must be specified. In our kinematic wave example, we need topography. We'll take elevation, and calculate slope. Therefore, elevation will be an input variable, and slope an output (though we won't modify it). A more advanced model might also need spatially variable precipitation, infiltration capacity, and/or roughness, but here we'll treat these as uniform.\n",
    "\n",
    "The `_input_var_names` should be a tuple of strings. These strings give the formal names for the fields, and should use standard Landlab semantics. In brief, Landlab's naming scheme follows the pattern of the CSDMS standard names, with the quantity at the end of the string, and any descriptors beforehand, separated from the quantity by two underscores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "_input_var_names = (\n",
    "    'topographic__elevation',\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Header information: `_output_var_names`\n",
    "*Output variables* are those arrays (fields) that the component solves for as part of its normal operation. In our kinematic wave example, the component will create fields for topographic gradient, water depth, velocity, and discharge per unit width (a.k.a., specific discharge).\n",
    "\n",
    "Like `_input_var_names`, `_output_var_names` should be either a tuple of strings with names based on <a href=\"https://github.com/landlab/landlab/wiki/Components#landlab-standard-naming-conventions\">standard Landlab semantics</a> ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "_output_var_names = (\n",
    "    'topographic__gradient', \n",
    "    'surface_water__depth',\n",
    "    'water__velocity',\n",
    "    'water__specific_discharge',\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Header information: `_var_units`\n",
    "This is a dictionary that defines the units of each variable (both input and output)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "_var_units = {\n",
    "    'topographic__elevation': 'm',\n",
    "    'topographic__gradient': 'm/m',\n",
    "    'surface_water__depth': 'm',\n",
    "    'water__velocity': 'm/s',\n",
    "    'water__specific_discharge': 'm2/s',\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Header information: `_var_mapping`\n",
    "This is a dictionary that describes the grid element to which each field is mapped. In our example, elevation and water depth (both scalars) are mapped to nodes, whereas velocity and discharge (both vectors) are mapped to links."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "_var_mapping = {\n",
    "    'topographic__elevation': 'node',\n",
    "    'topographic__gradient': 'link',\n",
    "    'surface_water__depth': 'node',\n",
    "    'water__velocity': 'link',\n",
    "    'water__specific_discharge': 'link',\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Header information: `_var_doc`\n",
    "This final piece of header data provides information about what each variable represents, encoded as a dictionary in which the keys are the field names and values are the descriptions, as strings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "_var_doc = {\n",
    "    'topographic__elevation':\n",
    "        'elevation of the ground surface relative to some datum',\n",
    "    'topographic__gradient':\n",
    "        'gradient of the ground surface',\n",
    "    'surface_water__depth':\n",
    "        'depth of water',\n",
    "    'water__velocity':\n",
    "        'flow velocity component in the direction of the link',\n",
    "    'water__specific_discharge':\n",
    "        'flow discharge component in the direction of the link',    \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Class with complete header information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from landlab import Component, FieldError\n",
    "import numpy as np\n",
    "\n",
    "class KinwaveOverlandFlowModel(Component):\n",
    "    \"\"\"\n",
    "    Calculate water flow over topography.\n",
    "    \n",
    "    Landlab component that implements a two-dimensional \n",
    "    kinematic wave model.\n",
    "    \n",
    "    Construction:\n",
    "    \n",
    "        KinwaveOverlandFlowModel(grid, [stuff to be added later])\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    grid : ModelGrid\n",
    "        A Landlab grid object.\n",
    "    [others to be added later]\n",
    "    \"\"\"\n",
    "\n",
    "    _name = 'KinwaveOverlandFlowModel'\n",
    "\n",
    "    _input_var_names = (\n",
    "        'topographic__elevation',\n",
    "    )\n",
    "\n",
    "    _output_var_names = (\n",
    "        'topographic__gradient',\n",
    "        'surface_water__depth',\n",
    "        'water__velocity',\n",
    "        'water__specific_discharge',\n",
    "    )\n",
    "\n",
    "    _var_units = {\n",
    "        'topographic__elevation': 'm',\n",
    "        'topographic__gradient': 'm/m',\n",
    "        'surface_water__depth': 'm',\n",
    "        'water__velocity': 'm/s',\n",
    "        'water__specific_discharge': 'm2/s',\n",
    "    }\n",
    "    \n",
    "    _var_mapping = {\n",
    "        'topographic__elevation': 'node',\n",
    "        'topographic__gradient': 'link',\n",
    "        'surface_water__depth': 'node',\n",
    "        'water__velocity': 'link',\n",
    "        'water__specific_discharge': 'link',\n",
    "    }\n",
    "\n",
    "    _var_doc = {\n",
    "        'topographic__elevation':\n",
    "            'elevation of the ground surface relative to some datum',\n",
    "        'topographic__gradient':\n",
    "            'gradient of the ground surface',\n",
    "        'surface_water__depth':\n",
    "            'depth of water',\n",
    "        'water__velocity':\n",
    "            'flow velocity component in the direction of the link',\n",
    "        'water__specific_discharge':\n",
    "            'flow discharge component in the direction of the link',    \n",
    "    }\n",
    "    def __init__(): # ignore this for now, we will add more stuff eventually. \n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The initialization method (`__init__`)\n",
    "Every Landlab component should have an `__init__` method. The parameter signature should start with a `ModelGrid` object as the first parameter. Following this are component-specific parameters. In our example, the parameters for the kinematic wave model include: precipiation rate, precipitation duration, infiltration rate, and roughness coefficient (Manning's n).\n",
    "\n",
    "Two things a component `__init__` method common does are (1) store the component's parameters as class attributes, and (2) create the necessary fields. When creating grid fields, it is important to first check to see whether a field with the same name (and mapping) already exists. For example, a driver or another component might have already created `topographic__elevation` when our kinematic wave component is initialized.\n",
    "\n",
    "(NOTE: IN LANDLAB 2.0, THE COMPONENT BASE CLASS PROVIDES A METHOD FOR AUTOMATICALLY CHECKING FOR REQUIRED INPUT FIELDS)\n",
    "\n",
    "There are a couple of different approaches to creating and working with fields. In the example below, we use the straightforward approach of creating them one at a time by name. We keep a reference to each field as a variable; for example, `self.elev` refers to the grid field `topographic_elevation`. It is also possible to automate the process by looping over all names in, say, `_var_mapping`, checking each to see whether it exists, and creating it if not. The `Component` method `initialize_output_fields()` may help with this, which will already be a member of your component because it inherits from `Component`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def __init__(self, grid, precip_rate=1.0, precip_duration=1.0, \n",
    "             infilt_rate=0.0, roughness=0.01):\n",
    "    \"\"\"Initialize the KinwaveOverlandFlowModel.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    grid : ModelGrid\n",
    "        Landlab ModelGrid object\n",
    "    precip_rate : float, optional (defaults to 1 mm/hr)\n",
    "        Precipitation rate, mm/hr\n",
    "    precip_duration : float, optional (defaults to 1 hour)\n",
    "        Duration of precipitation, hours\n",
    "    infilt_rate : float, optional (defaults to 0)\n",
    "        Maximum rate of infiltration, mm/hr\n",
    "    roughness : float, defaults to 0.01\n",
    "        Manning roughness coefficient, s/m^1/3\n",
    "    \"\"\"\n",
    "\n",
    "    # Store grid and parameters and do unit conversion\n",
    "    self._grid = grid\n",
    "    self.precip = precip_rate / 3600000.0 # convert to m/s\n",
    "    self.precip_duration = precip_duration * 3600.0  # convert to s\n",
    "    self.infilt = infilt_rate / 3600000.0 # convert to m/s\n",
    "    self.vel_coef = 1.0 / roughness  # do division now to save time\n",
    "\n",
    "    # Create fields...\n",
    "    #   Elevation\n",
    "    if 'topographic__elevation' in grid.at_node:\n",
    "        self.elev = grid.at_node['topographic__elevation']\n",
    "    else:\n",
    "        raise FieldError(\n",
    "            'A topography is required as a component input!')\n",
    "    #   Slope\n",
    "    if 'topographic__gradient' in grid.at_link:\n",
    "        self.slope = grid.at_link['topographic__gradient']\n",
    "    else:\n",
    "        self.slope = grid.add_zeros('link', 'topographic__gradient')\n",
    "        # this field is \"nice to have\", but we build it if necessary\n",
    "    #  Water depth\n",
    "    if 'surface_water__depth' in grid.at_node:\n",
    "        self.depth = grid.at_node['surface_water__depth']\n",
    "    else:\n",
    "        self.depth = grid.add_zeros('node', 'surface_water__depth')\n",
    "    #  Velocity\n",
    "    if 'water__velocity' in grid.at_link:\n",
    "        self.vel = grid.at_link['water__velocity']\n",
    "    else:\n",
    "        self.vel = grid.add_zeros('link', 'water__velocity')\n",
    "    #  Discharge\n",
    "    if 'water__specific_discharge' in grid.at_link:\n",
    "        self.disch = grid.at_link['water__specific_discharge']\n",
    "    else:\n",
    "        self.disch = grid.add_zeros('link',\n",
    "                                    'water__specific_discharge')\n",
    "\n",
    "    # Calculate the ground-surface slope (assume it won't change)\n",
    "    self.slope[self._grid.active_links] = \\\n",
    "        self._grid.calc_grad_at_link(self.elev)[self._grid.active_links]\n",
    "    self.sqrt_slope = np.sqrt( self.slope )\n",
    "    \n",
    "    # Store the current time (initializing it to zero)\n",
    "    self.current_time = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The \"go\" method, `run_one_step()`\n",
    "Every Landlab component will have a method that implements the component's action. The go method can have any name you like, but the preferred practice for time-advancing components is to use the standard name `run_one_step()`. Landlab assumes that if a component has a method with this name, it will (a) be the primary \"go\" method, and (b) will be fully standardized as described here.\n",
    "\n",
    "The `run_one_step` method should take either zero or one argument. If there is an argument, it should be a duration to run, `dt`; i.e., a timestep length. If the component does not evolve as time passes, this argument may be missing (see, e.g., the FlowRouter, which returns a steady state flow pattern independent of time).\n",
    "\n",
    "The first step in the algorithm in the example below is to calculate water depth *at the links*, where we will be calculating the water discharge. In this particular case, we'll use the depth at the upslope of the two nodes. The grid method to do this, `map_value_at_max_node_to_link`, is one of many mapping functions available.\n",
    "\n",
    "We then calculate velocity using the Manning equation, and specific discharge by multiplying velocity by depth.\n",
    "\n",
    "Mass balance for the cells around nodes is computed using the `calc_flux_div_at_node` grid method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "    def run_one_step(self, dt):\n",
    "        \"\"\"Calculate water flow for a time period `dt`.\n",
    "        \"\"\"\n",
    "        # Calculate water depth at links.\n",
    "        H_link = self.grid.map_value_at_max_node_to_link(\n",
    "                'topographic__elevation', 'surface_water__depth')\n",
    "\n",
    "        # Calculate velocity\n",
    "        self.vel = self.vel_coef * H_link**0.66667 * self.sqrt_slope\n",
    "\n",
    "        # Calculate discharge\n",
    "        self.disch = H_link * self.vel\n",
    "\n",
    "        # Flux divergence\n",
    "        dqda = self.grid.calc_flux_div_at_node(self.disch)\n",
    "\n",
    "        # Rate of change of water depth\n",
    "        if self.current_time < self.precip_duration:\n",
    "            ppt = self.precip\n",
    "        else:\n",
    "            ppt = 0.0\n",
    "        dHdt = ppt - self.infilt - dqda\n",
    "\n",
    "        # Update water depth\n",
    "        self.depth[self.grid.core_nodes] += dHdt[self.grid.core_nodes] * dt\n",
    "\n",
    "        # Somewhat crude numerical hack: prevent negative water depth\n",
    "        self.depth[np.where(self.depth < 0.0)[0]] = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changes to boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes, as in this example, it proves convenient to hard-code assumptions about boundary conditions into the `__init__` method. In this case, note that the final two lines of `__init__` invoke `grid.active_links` - so if the boundary conditions are updated after we initialize, it's possible `self.depth` will not be correct.\n",
    "\n",
    "We can resolve this issue by creating an additional component method that updates these components that can be called if the boundary conditions change. Whether the boundary conditions have changed can be assessed with a grid method called `bc_set_code`. This is an `int` which will change if the boundary conditions change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "    def __init__(self, grid, precip_rate=1.0, precip_duration=1.0, \n",
    "                 infilt_rate=0.0, roughness=0.01):\n",
    "        \"\"\"Initialize the KinwaveOverlandFlowModel.\n",
    "        ...\n",
    "        \"\"\"\n",
    "        # Store grid and parameters and do unit conversion\n",
    "        self._grid = grid\n",
    "        self._bc_set_code = self.grid.bc_set_code\n",
    "        # ...\n",
    "\n",
    "    def updated_boundary_conditions(self):\n",
    "        \"\"\"Call if boundary conditions are updated.\n",
    "        \"\"\"\n",
    "        # Calculate the ground-surface slope (assume it won't change)\n",
    "        self.slope[self.grid.active_links] = \\\n",
    "            self.grid.calc_grad_at_link(self.elev)[self.grid.active_links]\n",
    "        self.sqrt_slope = np.sqrt(self.slope)\n",
    "\n",
    "    def run_one_step(self, dt):\n",
    "        \"\"\"Calculate water flow for a time period `dt`.\n",
    "        \"\"\"\n",
    "        if self._bc_set_code != self.grid.bc_set_code:\n",
    "            self.updated_boundary_conditions()\n",
    "            self._bc_set_code = self.grid.bc_set_code\n",
    "        # Calculate water depth at links.\n",
    "        # ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The complete component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from landlab import Component, FieldError\n",
    "import numpy as np\n",
    "\n",
    "class KinwaveOverlandFlowModel(Component):\n",
    "    \"\"\"\n",
    "    Calculate water flow over topography.\n",
    "    \n",
    "    Landlab component that implements a two-dimensional \n",
    "    kinematic wave model.\n",
    "    \n",
    "    Construction:\n",
    "    \n",
    "        KinwaveOverlandFlowModel(grid, precip_rate=1.0, \n",
    "                                 precip_duration=1.0, \n",
    "                                 infilt_rate=0.0,\n",
    "                                 roughness=0.01, **kwds)\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    grid : ModelGrid\n",
    "        A Landlab grid object.\n",
    "    precip_rate : float, optional (defaults to 1 mm/hr)\n",
    "        Precipitation rate, mm/hr\n",
    "    precip_duration : float, optional (defaults to 1 hour)\n",
    "        Duration of precipitation, hours\n",
    "    infilt_rate : float, optional (defaults to 0)\n",
    "        Maximum rate of infiltration, mm/hr\n",
    "    roughnes : float, defaults to 0.01\n",
    "        Manning roughness coefficient, s/m^1/3\n",
    "        \n",
    "    Examples\n",
    "    --------\n",
    "    >>> from landlab import RasterModelGrid\n",
    "    >>> rg = RasterModelGrid((4, 5), 10.0)\n",
    "    >>> kw = KinwaveOverlandFlowModel(rg)\n",
    "    >>> kw.vel_coef\n",
    "    100.0\n",
    "    >>> rg.at_node['surface_water__depth']\n",
    "    array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
    "            0.,  0.,  0.,  0.,  0.,  0.,  0.])    \n",
    "    \"\"\"\n",
    "\n",
    "    _name = 'KinwaveOverlandFlowModel'\n",
    "\n",
    "    _input_var_names = (\n",
    "        'topographic__elevation',\n",
    "    )\n",
    "\n",
    "    _output_var_names = (\n",
    "        'topographic__gradient',\n",
    "        'surface_water__depth',\n",
    "        'water__velocity',\n",
    "        'water__specific_discharge',\n",
    "    )\n",
    "\n",
    "    _var_units = {\n",
    "        'topographic__elevation' : 'm',\n",
    "        'topographic__slope' : 'm/m',\n",
    "        'surface_water__depth' : 'm',\n",
    "        'water__velocity' : 'm/s',\n",
    "        'water__specific_discharge' : 'm2/s',\n",
    "    }\n",
    "\n",
    "    _var_mapping = {\n",
    "        'topographic__elevation' : 'node',\n",
    "        'topographic__gradient' : 'link',\n",
    "        'surface_water__depth' : 'node',\n",
    "        'water__velocity' : 'link',\n",
    "        'water__specific_discharge' : 'link',\n",
    "    }\n",
    "\n",
    "    _var_doc = {\n",
    "        'topographic__elevation':\n",
    "            'elevation of the ground surface relative to some datum',\n",
    "        'topographic__gradient':\n",
    "            'gradient of the ground surface',\n",
    "        'surface_water__depth':\n",
    "            'depth of water',\n",
    "        'water__velocity':\n",
    "            'flow velocity component in the direction of the link',\n",
    "        'water__specific_discharge':\n",
    "            'flow discharge component in the direction of the link',    \n",
    "    }\n",
    "    \n",
    "    _TWO_THIRDS = 2./3.\n",
    "\n",
    "    def __init__(self, grid, precip_rate=1.0, precip_duration=1.0, \n",
    "                 infilt_rate=0.0, roughness=0.01):\n",
    "        \"\"\"Initialize the KinwaveOverlandFlowModel.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        grid : ModelGrid\n",
    "            Landlab ModelGrid object\n",
    "        precip_rate : float, optional (defaults to 1 mm/hr)\n",
    "            Precipitation rate, mm/hr\n",
    "        precip_duration : float, optional (defaults to 1 hour)\n",
    "            Duration of precipitation, hours\n",
    "        infilt_rate : float, optional (defaults to 0)\n",
    "            Maximum rate of infiltration, mm/hr\n",
    "        roughnes : float, defaults to 0.01\n",
    "            Manning roughness coefficient, s/m^1/3\n",
    "        \"\"\"\n",
    "\n",
    "        # Store grid and parameters and do unit conversion\n",
    "        self._grid = grid\n",
    "        self._bc_set_code = self.grid.bc_set_code\n",
    "        self.precip = precip_rate / 3600000.0 # convert to m/s\n",
    "        self.precip_duration = precip_duration * 3600.0  # h->s\n",
    "        self.infilt = infilt_rate / 3600000.0 # convert to m/s\n",
    "        self.vel_coef = 1.0 / roughness  # do division now to save time\n",
    "\n",
    "        # Create fields...\n",
    "        #   Elevation\n",
    "        if 'topographic__elevation' in grid.at_node:\n",
    "            self.elev = grid.at_node['topographic__elevation']\n",
    "        else:\n",
    "            raise FieldError(\n",
    "                'A topography is required as a component input!')\n",
    "        #  Water depth\n",
    "        if 'surface_water__depth' in grid.at_node:\n",
    "            self.depth = grid.at_node['surface_water__depth']\n",
    "        else:\n",
    "            self.depth = grid.add_zeros('node', 'surface_water__depth')\n",
    "        #   Slope\n",
    "        if 'topographic__gradient' in grid.at_link:\n",
    "            self.slope = grid.at_link['topographic__gradient']\n",
    "        else:\n",
    "            self.slope = grid.add_zeros('link', 'topographic__gradient')\n",
    "        #  Velocity\n",
    "        if 'water__velocity' in grid.at_link:\n",
    "            self.vel = grid.at_link['water__velocity']\n",
    "        else:\n",
    "            self.vel = grid.add_zeros('link', 'water__velocity')\n",
    "        #  Discharge\n",
    "        if 'water__specific_discharge' in grid.at_link:\n",
    "            self.disch = grid.at_link['water__specific_discharge']\n",
    "        else:\n",
    "            self.disch = grid.add_zeros('link',\n",
    "                                        'water__specific_discharge')\n",
    "\n",
    "        # Calculate the ground-surface slope\n",
    "        self.slope[self._grid.active_links] = \\\n",
    "            self._grid.calc_grad_at_link(self.elev)[self._grid.active_links]\n",
    "        self.sqrt_slope = np.sqrt( self.slope )\n",
    "        self.sign_slope = np.sign( self.slope )\n",
    "\n",
    "        # Store the current time (initializing it to zero)\n",
    "        self.current_time = 0.0\n",
    "        \n",
    "    def updated_boundary_conditions(self):\n",
    "        \"\"\"Call if boundary conditions are updated.\n",
    "        \"\"\"\n",
    "        self.slope[self.grid.active_links] = \\\n",
    "            self.grid.calc_grad_at_link(self.elev)[self.grid.active_links]\n",
    "        self.sqrt_slope = np.sqrt(self.slope)\n",
    "        self.sign_slope = np.sign( self.slope )\n",
    "\n",
    "    def run_one_step(self, dt):\n",
    "        \"\"\"Calculate water flow for a time period `dt`.\n",
    "        \"\"\"\n",
    "\n",
    "        if self._bc_set_code != self.grid.bc_set_code:\n",
    "            self.updated_boundary_conditions()\n",
    "            self._bc_set_code = self.grid.bc_set_code\n",
    "\n",
    "        # Calculate water depth at links\n",
    "        H_link = self._grid.map_value_at_max_node_to_link(\n",
    "                'topographic__elevation', 'surface_water__depth')\n",
    "\n",
    "        # Calculate velocity\n",
    "        self.vel = -self.sign_slope * self.vel_coef * H_link**0.666667 \\\n",
    "                    * self.sqrt_slope\n",
    "\n",
    "        # Calculate discharge\n",
    "        self.disch = H_link * self.vel\n",
    "\n",
    "        # Flux divergence\n",
    "        dqda = self._grid.calc_flux_div_at_node(self.disch)\n",
    "\n",
    "        # Rate of change of water depth\n",
    "        if self.current_time < self.precip_duration:\n",
    "            ppt = self.precip\n",
    "        else:\n",
    "            ppt = 0.0\n",
    "        dHdt = ppt - self.infilt - dqda\n",
    "\n",
    "        # Update water depth\n",
    "        self.depth[self._grid.core_nodes] += dHdt[self._grid.core_nodes] * dt\n",
    "\n",
    "        # Somewhat crude numerical hack: prevent negative water depth\n",
    "        self.depth[np.where(self.depth < 0.0)[0]] = 0.0\n",
    "        \n",
    "        self.current_time += dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The discharge from node 6 to node 5 should be -0.000278 m2/s:\n",
      "-0.000277754737589\n",
      "The discharge from node 5 to node 4 should be -0.000556 m2/s:\n",
      "-0.000555461149416\n"
     ]
    }
   ],
   "source": [
    "from landlab import RasterModelGrid\n",
    "nr = 3\n",
    "nc = 4\n",
    "rg = RasterModelGrid((nr, nc), 10.0)\n",
    "rg.add_empty('node', 'topographic__elevation')\n",
    "rg.at_node['topographic__elevation'][:] = rg.x_of_node.copy()\n",
    "kinflow = KinwaveOverlandFlowModel(rg, precip_rate=100.0, precip_duration=100.0)\n",
    "for i in range(100):\n",
    "    kinflow.run_one_step(1.0)\n",
    "print('The discharge from node 6 to node 5 should be -0.000278 m2/s:')\n",
    "print(kinflow.disch[8])\n",
    "print('The discharge from node 5 to node 4 should be -0.000556 m2/s:')\n",
    "print(kinflow.disch[7])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll test the component on a larger grid and a larger domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "nr = 62\n",
    "nc = 42\n",
    "rg = RasterModelGrid((nr, nc), 10.0)\n",
    "rg.add_empty('node', 'topographic__elevation')\n",
    "rg.at_node['topographic__elevation'] = 0.01 * rg.y_of_node\n",
    "kinflow = KinwaveOverlandFlowModel(rg, precip_rate=100.0, precip_duration=100.0)\n",
    "for i in range(1800):\n",
    "    kinflow.run_one_step(1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the topography:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPIAAAEHCAYAAACUQBu6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAVUElEQVR4nO3df6xk5V3H8fenUKCthYVuwc3u6rZ2o22alJKbFlN/YFGzYOPyBxj8iXXj/iEY/F1qY4zGGBpNaZso8QZqF61uEW3YNKSVbEsak4JdSsXCqmwR2euubK8FiiDs3p2Pf8xzy7g7c+/M7Dn3zDn7eZGTmfPMM2efe8J3nuc85znPI9tERLu9oukCRMSpSyBHdEACOaIDEsgRHZBAjuiABHJEB5zZdAFOxfr1671ly5amixEd8MQTT7C4uKhTOca2bdu8uLg4Vt4HH3zws7a3rZRH0jrgNuCtgIFfsP3FYXlbHchbtmxh3759TRcjOmBubu6Uj7G4uDj2/4+S1o+R7SPAZ2xfLeks4NWjMtbatJa0TtJdkv5F0n5J3yvpAkn3SnqsvJ5f8krSRyUdkPSwpEvqLFtEPTzmtjJJ5wI/ANwOYPuo7WdG5a/7Gnn5F+V7gLcB+4GbgL22twJ7yz7AFcDWsu0Ebq25bBGVs3tjbWN4I/B14M8lPSTpNkmvGZW5tqb1wC/Kz0P/FwU4Kmk7cFnJtgu4D3gfsB24w/0xo/eX2nyD7cOj/xVjL9X0F8TppYqhyuPVtsV6SYPt8Hnb8wP7ZwKXAL9s+wFJH6Ff6f3OsIPVeY08+IvyNuBB4EbgouXgtH1Y0oUl/0bg4MD3F0raCoEcMVsmeHZh0fZKF+YLwILtB8r+Xbzcej1JnU3r5V+UW22/HXh+pYIAw3oMTzorknZK2idp39e/Pl4PYcTaqeYa2fZ/AQclfXdJuhx4dFT+OgN52C/KJcBTkjYAlNcjA/k3D3x/E3DoxIPanrc9Z3vu9a8fp+MvYhxVPQVYTSAXvwx8QtLDwMXAH47KWFvT2vZ/SToo6btt/ysv/6I8ClwH3Fxe7y5f2QPcIGk38E7g2ZWvj/t6vVwjx+yo8rFg218BxrovVvd95OVflLOAx4H30m8F3ClpB/AkcE3Jew9wJXAAeKHkjWgRA2P1SFeu1kBe4Rfl8iF5DVxfZ3ki6mRXWyNPotUju+zcfopqVBd/CeSIDkggR7Sc07SOaFZVAdjBzq765Ro5Zktq5IjWm2iwR6USyBGVStN6CqbXO9Z0IaITqqlJ07SO6IQEckTLedxJAyqXQI6gyiZxauQpGDvXyDFLEsgRrZaHJiI6IfeRp2Ln9lNUparbT8crOc6kWh3IEbMnNXJEy+Xpp4iG5fZTg/L0U8yaBHJE66VpHdF6BtJrPYXcfoqq5OmniPiWBHJEBySQIxpTRZO4P896AnlimaA+Zk91gSzpCeA5+j1oSystw1prIA8riKQLgE8CW4AngJ+w/bQkAR+hv/7TC8DP2/5yneWLqFoNY61/yPaq6wfXuazqYEEuHvg1uQnYa3srsJeX10y+Athatp3ArWtQtogKjbukavXN7yaa1tuBy8r7XcB9wPtK+h1lMbf7Ja2TtGHlpVVNr/dSrYWN08VMDtE08PeSDPyZ7flRGesO5GEFuWg5OG0flnRhybsRODjw3YWStuoayRGzYoI5u9ZL2jewPz8kUN9l+1CJkXsl/YvtLww7WN2BfFJBVsirIWkn/bxJ2km/6c2mTd9eTSkjKjN2jby4UucVgO1D5fWIpE8B7wCGBnKt18iDBQGWC/KUpA0A5fVIyb4AbB74+ibg0JBjztuesz23fv35dRY/TiPV3Dbqz6I5zrYaSa+R9Nrl98CPAl8dlb+2Grn846+w/dxAQX4f2ANcB9xcXu8uX9kD3CBpN/BO4NmVr48hk+/F7KlsOtyLgE/1b+ZwJvBXtj8zKnOdTeuhBZH0JeBOSTuAJ4FrSv576N96OkD/9tN7ayxbROWqnHzP9uPA28bNX1sgjyqI7f8GLh+SbuD6usoTUb9MvjeVTL4X1ZnJ209ja3UgR8yaLBkT0Xomy6pGNKqqTqoE8hRy+ylmTa6RIzoggRzRclkfeUq5/RRVye2niAD6NXKmw43ogNTIEY2pZvK9zGs9Fds418gxSxoK5LWYsysiatbqGjli1riXpvXkbHpLzdy3i46pbIaQBHJEuzX3OHICOaJSqZEj2q+hOG5/IPt4rpFjhqRGjmi/dHZFtF06u6Zjk9tPUYnKatLUyBHtZtLZFdEBTo0c0ah2t6xbHsg2zjVyzAoDFY61lnQGsA/4T9vvWSlv7U8/STpD0kOSPl323yDpAUmPSfqkpLNK+tll/0D5fEvdZYuomsf8b0w3AvvHybgWjzGeWJgPArfY3go8Dewo6TuAp22/Cbil5ItoF4+5rULSJuDHgNvG+WdrDeQTC6P+0ozvBu4qWXYBV5X328s+5fPLS/6I+lU2957H21b3YeC3GHPpirqvkZcL89qy/zrgGdtLZX8B2FjebwQOAthekvRsyb848uiG3lIzk51FDDNBZ9d6SfsG9udtzwNIeg9wxPaDki4b52B1LnQ+rDDDaliP8dngcXcCOwE2blhfQUkjKuKJJhZYtD034rN3AT8u6UrgHOBcSX9p+2dGHazOpvVyYZ4AdtNvUn8YWCdp+QdkE3CovF8ANgOUz88DvnHiQW3P256zPXfB+efWWPw4nUzQAbXiUapoWtt+v+1NtrcA1wKfWymIod6Fzt8PvB+g1Mi/YfunJf0NcDX94L4OuLt8ZU/Z/2L5/HNeZdycyQwhMWNOo/vI7wN2S/oD4CHg9pJ+O/AXkg7Qr4mvbaBsEVPrD9GsNpJt3wfct1q+NQnkwcLYfhx4x5A8LwLXrEV5ImqRp58iGlbZEM2MtZ6cyRDNmC2ZDjei/fLQRESTqojABh9IbncgZ4L6mDWpkSPaLitNRHRC1n6KaFBFSz+laT2V3H6KWZOmdUS7ZRbNiKZlXuvm2KZ3LE3rmB2pkSPazs4QzYguyH3kiCZVOfleA9odyJl8L2ZMrpEj2i4PTUQ0q7plVas5zKRaHciZfC9mTcZaR7Sc8/RTRAfkoYmIjkiNPIU8/RQzJk3riA5IZ1dEk6qafK+hBuLIQJZ0D/BLtp9Yu+JMyOTpp5gtDTWtV1qN8ePA30v6gKRXTnpgSedI+kdJ/yTpEUm/V9LfIOkBSY9J+qSks0r62WX/QPl8yxR/T0SjqlrnfFT8jDIykG3fCbwdOBfYJ+k3JP3a8jbG3/QS8G7bbwMuBrZJuhT4IHCL7a3A08COkn8H8LTtNwG3lHwRLVLNsqrFqPgZarX1kY8BzwNnA689YVv5T+r7n7L7yrKZ/jrJd5X0XcBV5f32sk/5/HJJwxY/j6hcVZPvVRXHK8TPUCtdI28DPkR/3eJLbL8wyd9UjnEG8CDwJuBPgK8Bz9heKlkWgI3l/UbgYPkjliQ9C7wOWBx1fNv0jubpp5ghver6bE6MH9sPjMq7Uq/1B4BrbD8ybUFsHwculrQO+BTw5mHZyuuw2vekXyBJO4GdABvWnz9t0SIqZ8Djx/F6SfsG9udtz/+/450QP5Leavurww42MpBtf//YRVqF7Wck3QdcCqyTdGaplTcBh0q2BWAzsCDpTOA8+guen3iseWAe4C3f9R0NDYiLGGKyxxgXbc+NddiX42cbMFkgnypJrweOlUK8Cvhh+h1YnweuBnYD1wF3l6/sKftfLJ9/zqsNk7FxJhaIKlR026iqu08rxM9QdQ4I2QDsKu38VwB32v60pEeB3ZL+AHgIuL3kvx34C0kH6NfE19ZYtoh6VHcfeWj8jMpcWyDbfpj+7asT0x8H3jEk/UXgmrrKE1G/sW8trX6kEfEzSoZoRkA1QysNPp6x1pPLEM2YMZl8L6IL8hhjRMs5zyNPxWlaR0UqC8BZe4wxIibTH9mVGjmi3bKIW0TDKlseOYE8ORsfyxDNmCG5/RTRAamRI1rO6eyaTs8cfylN66hAVWOkE8gRLTeL0+FGxKSyiFtEoyqLvzStJ2eDM/lezJLcfopoN+ehiYhuyMQC07DpHc3TT1GByhZxSyBHtFx6rSO6IfeRI5pTVUWaGnkahl6efopZYSCdXRHtlxo5ouX6Sz8lkCdmm2PH07SOClQUgE3dDF1tofOpSdos6fOS9kt6RNKNJf0CSfdKeqy8nl/SJemjkg5IeljSJXWVLaIW7t9+GmerWm2BDCwBv277zfSXU71e0luAm4C9trcCe8s+wBXA1rLtBG6tsWwRtagqkEdVhKPUFsi2D9v+cnn/HLAf2AhsB3aVbLuAq8r77cAd7ruf/jrKG+oqX8SgKmrJ/sAuj7WNYVRFONSaXCNL2kJ/ZbkHgItsH4Z+sEu6sGTbCBwc+NpCSTs86riGXCPHTKmq2VxiZDlOnpO0XBE+Oix/7YEs6duAvwV+xfY3JY3MOiTtpLMiaSf9pjcXnntuVcWMqEQd178nVIRD1XmNjKRX0g/iT9j+u5L81HKTubweKekLwOaBr28CDp14TNvztudsz5336lfXV/iIKXjMDVgvad/AtnPY8U6sCEf9u7XVyOpXvbcD+21/aOCjPcB1wM3l9e6B9Bsk7QbeCTy73AQfxTZLaVpHBaqoRyfskV60PbdShhEV4VB1Nq3fBfws8M+SvlLSfpt+AN8paQfwJHBN+ewe4ErgAPAC8N4ayxZRizE7sla1QkU4VG2BbPsfGH7dC3D5kPwGrq+rPBFrocJr5KEVoe17hmVu9ciuiKpU2Ntc1XFWqghP0upAzhDNmCUZax3REU2NtU4gR1SlpnHU42h1IBty+ykqUcntJ6DXa6ZObnUgR8yahuanTyBHVClN64gGzdrtp0m1OpBz+ylmicd/RLFyrQ7kiFmTa+SIDkiv9RRsc3RpqeliRAdU1STONXJEyy1P9dOEBHJEVTKyK6JZVYVfOrumYOBYrpFjRhg4ns6uiPZL0zqiAxLIU8jIrqhKJRPUp7MrohsysUBEB6RGjmhQVWs/pdd6Cr0M0YwZkxo5ou3S2RXRfhlrPSXbGdkVlWj7DCG1rcYo6WOSjkj66kDaBZLulfRYeT2/pEvSRyUdkPSwpEvqKldEnZbvJa+2Va3OZVU/Dmw7Ie0mYK/trcDesg9wBbC1bDuBW2ssV0QtbHO81xtrW82winAltQWy7S8A3zgheTuwq7zfBVw1kH6H++4H1i2voRyxFqqqJXtl3q7VtjF8nJMrwpHW+hr5ouU1j20flnRhSd8IHBzIt1DSVl0f+cVjx2opaMQ0KrzW/oKkLePmn5XOrmGrzg09I2Vl950A573qVXWWKWIiTfZa13mNPMxTy03m8nqkpC8AmwfybQIODTuA7Xnbc7bnXnP22bUWNmJSE3R2rZe0b2DbeSr/7lrXyHuA64Cby+vdA+k3SNoNvBN4drkJvpJMvhdVqaRJXDq7xrRoe+7U/9G+2gJZ0l8Dl9H/5VkAfpd+AN8paQfwJHBNyX4PcCVwAHgBeG9d5YqoSycHhNj+yREfXT4kr4Hr6ypLxFqpqrNrWEVo+/ZR+WelsyuiUZXUpBUuGbNCRThUqwO5B7lGjplh8vRTRCd07ho54nSTiQWmZJuXMrIrKlBJPZplVSO6IYEc0XIGnKZ1RHOqqklTI08h18gxU3KNHNF+BpYaWvkkgRxRoQwImUKv1+OFl15quhjRAb0KOqmcpnVEN1TxgzCNBHJERTr5GGNEm1Q1sUA6u6bQs/nfo0ebLkYEUMZap0aOaL9cI0e0XHqtp9Szc/spKlFVkzg1ckTL5XnkiA6wzbH0Wkc0p6omcZrWUzje6/H8iy82XYwI4OXVGJvQ6kCOmCW5Ro7ogtTI0+n1enzz+eebLkZ0QBUBaHKNDICkbcBHgDOA22zf3HCRIsZmm5caGjI8M4Es6QzgT4Afob/M6pck7bH9aLMlixiPbZYaqpHXen3klbwDOGD7cdtHgd3A9obLFKeJ4xXN/Xb8+PGxtnFI2ibpXyUdkHTTSnlnpkYGNgIHB/YX6K+VPNLSsWMsHhq6HnrEmnOvx7GKhgxP2kKdpUDWkLSTBsCWld13Qv9COmJW2OZYddfI32qhAkhabqHOfCAvAJsH9jcBJ1W3tueBeYCzpWYeNYkYotfr8cJzz1V1uIlaqLMUyF8Ctkp6A/CfwLXAT630haOw+CT8B7AeWKy/iK2T8zLcsPPynad60KPw2Sf7xx7HOZL2DezPl0pq2Vgt1GUzE8i2lyTdAHyWfqv5Y7YfWeU7rweQtM/23BoUs1VyXoar67zY3lbh4cZqoS6bmUAGsH0PcE/T5YiYARO1UGcqkCOib9IWalcCeX71LKelnJfhWnFeJmmhqqklLiKiOrM0sisiptTqQJ5kCFvXSNos6fOS9kt6RNKNJf0CSfdKeqy8nl/SJemj5Vw9LOmSZv+Cekk6Q9JDkj5d9t8g6YFyXj4p6aySfnbZP1A+39JkuafV2kAeGMJ2BfAW4CclvaXZUq2pJeDXbb8ZuBS4vvz9NwF7bW8F9pZ96J+nrWXbCdy69kVeUzcC+wf2PwjcUs7L08COkr4DeNr2m4BbSr7WaW0gc5o/ZGH7sO0vl/fP0f+fdiP9c7CrZNsFXFXebwfucN/9wDpJG9a42GtC0ibgx4Dbyr6AdwN3lSwnnpfl83UXcHnJ3yptDuRhQ9g2NlSWRpXm4NuBB4CLbB+GfrADF5Zsp9P5+jDwW8DyM4WvA56xvVT2B//2b52X8vmzJX+rtDmQJxrC1lWSvg34W+BXbH9zpaxD0jp3viS9Bzhi+8HB5CFZPcZnrdHm+8gTDWHrIkmvpB/En7D9dyX5KUkbbB8uTecjJf10OV/vAn5c0pXAOcC59GvodZLOLLXu4N++fF4WJJ0JnAd8Y+2LfWraXCN/awhb6YG8FtjTcJnWTLmOux3Yb/tDAx/tAa4r768D7h5I/7nSe30p8OxyE7xLbL/f9ibbW+j/P/E52z8NfB64umQ78bwsn6+rS/7W1cjYbu0GXAn8G/A14ANNl2eN//bvo98EfBj4StmupH99txd4rLxeUPKLfi//14B/Buaa/hvW4BxdBny6vH8j8I/AAeBvgLNL+jll/0D5/I1Nl3uaLSO7IjqgzU3riCgSyBEdkECO6IAEckQHJJAjOiCB3FHl6ah/l3RB2T+/7J/yJHMxexLIHWX7IP0nnJbXz7qZ/kyN/9FcqaIuuY/cYWUI54PAx4BfBN7u/pNi0TFtHmsdq7B9TNJvAp8BfjRB3F1pWnffFcBh4K1NFyTqk0DuMEkX018E7FLgV7s6kUAkkDurPB11K/3nlJ8E/gj442ZLFXVJIHfXLwJP2r637P8p8D2SfrDBMkVN0msd0QGpkSM6IIEc0QEJ5IgOSCBHdEACOaIDEsgRHZBAjuiABHJEB/wfnIa/cFyY6OgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from landlab.plot import imshow_grid\n",
    "imshow_grid(rg, 'topographic__elevation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The steady solution should be as follows. The unit discharge at the bottom edge should equal the precipitation rate, 100 mm/hr, times the slope length.\n",
    "\n",
    "The slope length is the distance from the bottom edge of the bottom-most row of cells, to the top edge of the top-most row of cells. The base row of nodes are at y = 0, and the cell edges start half a cell width up from that, so y = 5 m. The top of the upper-most row of cells is half a cell width below the top grid edge, which is 610 m, so the top of the cells is 605 m. Hence the interior (cell) portion of the grid is 600 m long.\n",
    "\n",
    "Hence, discharge out the bottom should be 100 mm/hr x 600 m = 0.1 m/hr x 600 m = 60 m2/hr. Let's convert this to m2/s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.016666666666666666"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_out = 0.1 * 600 / 3600.0\n",
    "q_out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The water depth should be just sufficient to carry this discharge with the given slope and roughness. We get this by inverting the Manning equation:\n",
    "\n",
    "$$q = (1/n) H^{5/3} S^{1/2}$$\n",
    "\n",
    "$$H^{5/3} = n q S^{-1/2}$$\n",
    "\n",
    "$$H = (n q)^{3/5} S^{-3/10}$$\n",
    "\n",
    "The slope gradient is 0.01 (because we set elevation to be 0.01 times the y coordinate). The discharge, as we've already established, is about 0.0167 m2/s, and the roughness is 0.01 (the default value). Therefore,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.021559063053142917"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 0.01\n",
    "q = 0.0167\n",
    "S = 0.01\n",
    "H_out = (n * q)**0.6 * S**-0.3\n",
    "H_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEGCAYAAAB7IBD2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAb5ElEQVR4nO3df/BddX3n8ecrPwCLshESGCo/1Wil7vAjGWGH7q6FxU2oBTrWGrSwMLZpO7DColtFx112R2ewa4Gy1MxEpDWO8mNRxtSiGCmMtSOUJLAgBgdEwGAWEuRXQZJ8w3v/OJ+bnO/lnu899+ac7z333NfDOfP93nPP/Xw/3++YN5/fb0UEZmaDmjPqCpjZeHLwMLOhOHiY2VAcPMxsKA4eZjaUeaOuwN5YuHBhHHnkUaOuhrXA448/xrZt27Q3Zcw94MiIqV+VejZ+tfW2iFi2Nz9v1MY6eBx55FH8093rR10Na4GTT1y612XE1K/Y9+1/UOrZV+7764V7/QNHrNZui6QFkm6W9JCkTZL+jaQDJa2T9HD6+sb0rCRdLekRSfdLOqHOuplVT6A55a4WqPu3+CvgOxHxG8CxwCbgE8DtEbEYuD29BlgOLE7XSmBVzXUzq5aAOXPLXS1QW7dF0gHAvwPOA4iIHcAOSWcC706PfRm4E/g4cCawJrIlr3elVsuhEbGlzM977qWdldbf2m3B/vPrKVh7NWwyVupsebwZ2Ar8jaR7JV0raX/gkE5ASF8PTs+/Cfh57vOb0z2zMeFuS1XmAScAqyLieOAl9nRReukVsl+z8UbSSknrJa3fum1rNTU1q4pU7mqBOoPHZmBzRNydXt9MFkyeknQoQPr6dO75w3OfPwz4RXehEbE6IpZGxNJFCxfVVnmzgYmJannUNuYREf9P0s8lvT0ifgKcCvw4Xf8JuDx9/Wb6yFrgQkk3ACcCz5cd7wB4eceuSutv7VbPmEd7WhVl1L3O4z8DX5W0D/AocD5Za+cmSR8GngDen569FTgdeAR4OT1rNl5aMpNSRq3BIyLuA3qtvjm1x7MBXFBnfczqpdZ0ScoY6xWmedun3G2xERPutpjZkCao5TE5v6lZ7apd5yFpmaSfpC0br1nmIGlfSTem9++WdFS6f5qkDZIeSF9PyX1mSbr/SNoOonS/57aRmTh4mFVFwNy55a5+RUlzgb8m27ZxDHC2pGO6Hvsw8GxEvBW4Evhcur8N+N2I+NdkM5pfyX1mFdn2j85WkM7O3qJtI4Va023ZMfXqqKtgVuWYx7uARyLi0axY3UC2hePHuWfOBC5L398MXCNJEXFv7pkHgf0k7QscCBwQET9MZa4BzgK+ncp6d/pMfttIIbc8zCozULdlYWeldLpWdhVWZrvG7mciYgp4Hjio65n3AfdGxPb0/OaCMou2jRRqTcvDrBHKtzy2RcRMh4iU2a4x4zOSfpOsK/OeAcosrTXBY+eU889YA1Q321Jmu0bnmc2S5gH/CvglgKTDgFuAcyPip7nnDyso86nOLvaubSOF3G0xq0rZTXHlWif3AIslHZ1WaK8g28KRt5ZsQBTg94F/iIiQtAD4e+DSiPinzsOpO/KipJPSLMu5TN8e0ikrv22kkIOHWZUqOgwojWFcCNxGdojWTRHxoKT/KemM9NiXgIMkPQJcwp4ZkguBtwKflnRfujpjGH8GXEu2DeSnZIOlkO01O03Sw8Bp6fWMWtNtMRu9apenR8StZHu+8vf+W+77V9izNyz/zGeAzxSUuR54Z4/7z9Bj28hMWhM8du7yVK01gJenm9nAOud5TAgHD7PKeFftWHpll3fVWgP4PA8zG4rHPMxsYHK3xcyG5ZbH+NnuXbXWAHLwMLNBZacQOniY2aAkNMfBY+x4qtaawC0PMxuKg4eZDcXBw8wGJ3qf1dVSrQkeO171VK2NlpBbHmY2nDlzJmeFaa2/qaTHUoKZ+yStT/d6JpdR5uqUjOZ+SSfUWTezOkgqdZUsa9ikTwdJukPSv0i6Jvf8G3Ini90naZukq9J750namnvvj/rVbzZaHr8dEdtyrzvJZS5Pf5BPkOWHWM6eRDQnkiWnObHsD9nhw4Bs1Coc88glfTqN7ODieyStjYh83pbdSZ8krSA7Kf0DwCvAp8lODNt9alhEvAgcl/sZG4Bv5Mq7MSIuLFvHUbSxziRLKkP6elbu/prI3AUsSKc4m42NClseu5M+RcQOoJP0KS//b+lm4NSU9OmliPgBWRApqudistws/zjo79hRd/AI4LspX2YnqU1RcpkySW6QtLKTKGfrtq01Vt1sMJ0B05LBY7aSPhU5m6ylkc9Z8r40ZHCzpMOLPthRd7fl5Ij4RTq5eZ2kh2Z4tlRCmohYDawGWLJkqZO1WKMMsDy99qRPfawAzsm9/jvg+ojYLulPyVo0p/T8ZFJr8IiIX6SvT0u6hawpVpRcpkySm0LbPeZho6ZKF4ntVdKnGaspHQvMi4gNnXvp9PSOL7InaXah2rotkvaX9IbO92Qp735EcXKZtcC5adblJOD5TvfGbFxUOOYxdNKnEmWfDVzfVe/8+OIZZLliZlRny+MQ4Jb0h5oHfC0iviPpHuAmSR8GnmBP3olbgdPJktG8DJxfY93MalFVyyMipiR1kj7NBa7rJH0C1kfEWrKkT19RlvTpl2QBplOPx4ADgH0knQW8JzdT8wdk/9byPqIsmdRUKuu8fnWsLXhExKPAsT3u90wukyLmBcP+vFd2uttio1X1CtNhkz6l946aodw397h3KXDpIPXzClOzKk3O6nQHD7PKaLKWpzt4mFWoym5L07UmeHiq1hphcmJHe4KHWRO45WFmAxtkx2wbtCZ4vLLTK9Vt9Bw8zGwoTr1gZkNxy8PMBlftxrjGa03w8EliNmpZuslR12L2tCZ4mI2eZ1vMbEhzPGA6fl6Z8lStjZjcbTGzIQi3PMxsSJPU8pic/cNms6DJSZ/Se3emMjvJnQ6eqayZtKblsWPKU7U2YhWOedSR9CnnQxGxvuteUVmF3PIwq4gQc+bMKXWVUGvSpx56ljXTBxw8zCoklbsYbdKnv0ldlk/nAsTAZbWm27J9l6dqbfQGWCQ2qqRPH4qIJ1NalK+TJX5aM0xZbnmYVaVkq6NkfBkk6RNlkz5FxJPp64vA18i6R0OV5eBhVpFsb0tzkz5JmidpYfp+PvBeskRsA5cFLeq2mDVBVbMtdSR9Ah4HbkuBYy7wPbLUksxUVpHWBI+dnqq1BqhyhWlNSZ+WFDxfWFaR1gQPs5GbsPM8ah/zkDRX0r2SvpVeH51WsD2cVrTtk+4PvMLNrEk653lUNGDaeLMxYHoR0zNufw64MiIWA8+SrWyD3Ao34Mr0nNkYKTdY2pbWSa3dFkmHAb8DfBa4JC1IOQX4YHrky8BlwCqyFW6Xpfs3A9ek1XKlFnDs9Eli1gAtiQul1N3yuAr4c6DzL/sg4Lm0gg2mr5ortcJN0srOqryt27bWWXezwSgbMC1ztUFtwUPSe4GnI2JD/naPR6PEe3tuRKyOiKURsXTRwkUV1NSsGhWv82i8OrstJwNnSDod2I9szvkqYIGkeal1kV8111nhtrnsCrc8T9VaE7QlMJRRW8sjIi6NiMPSfPMKshVrHwLuIFvBBtmKtm+m7wde4WbWNJ5tqdfHyQZPHyEb0/hSuv8l4KB0/xLgNYefmDWduy0Vi4g7gTvT94+yZzNO/pmBV7iZNUqLWhVltGaFqadqbdSyw4AmJ3q0JniYNcGcCWp6OHiYVWiCYkd7gseUuy02YpqwjXGtCR5mTTBBQx4OHmZVmqQBUx9DaFYRkc24lPlfqfIqTvok6dck/b2khyQ9KOny3HvnSdqqPcmg/qhf/VrT8vDydGuCqhoeNSZ9+nxE3JHO0bld0vKI+HZ678aIuLBsHd3yMKtKydWlJQdVK0/6FBEvR8Qd6fsdwEay/WVDcfAwq1CFe1vqTPqEpAXA7wK3526/T9L9km6WdHjBR3drT7fFU7U2YmKgRWILJeXzxa6OiNVdxXWrIulTJy/L9cDVabsIwN8B10fEdkl/StaiOWWmcloTPMyaYIDZln4Z4wZJ+jToMRargYcj4qrOjYh4Jvf+FylxDKi7LWYVKdtlKdk4qTzpU1ZHfYYsyFzcdf/Q3MszmH7ucE9ueZhVqKq9LTUlfXoB+BTwELAxDdxeExHXAh+RdAYwlco6r18dWxM8ppzo2hqgyiViNSV96lnFiLgUuHSQ+rUmeJg1gfe2mNnAstmWUddi9rQmeHiFqY2cfBiQmQ3J3RYzG5i7LWY2NLc8xtAuL0+3Bpic0NGi4GE2ahLMnaB+i4OHWYXcbRlD7rZYE0xQ7CjeGCfp1s6xZmbWnxBzVO5qg5l21f4t8F1Jn5I0f9CCJe0n6Z8l/d90XuL/SPePTuctPpzOX9wn3e95HqPZ2Kh2V23jFQaPiLgJOJ5sZ956SR+TdEnnKlH2duCUiDgWOA5YJukksnMCroyIxcCzZOcwQu48RuBKSpwnYNY0TnS9x07gJWBf4A1A6YGFdK7Av6SX89MVZKcTfTDd/zJwGbCK7DzGy9L9m4Fr0nmMpbbLTnl5uo2YgLktCQxlFAYPScuAK8gOHDkhIl4etPB0AvQG4K1kJ0H/FHgunbcI089lnHYeo6TOeYzbuspcCawEOPyIIwatklmtJmimdsaWx6eA90fEg8MWHhG7gOPSYau3AO/o9Vj6Wuo8xnTO42qAJUuW+hAPaxQHDyAi/m1VPyQinpN0J3ASsEDSvNT6yJ/LOOx5jADs8mFANmKTlqu2tjNMJS1KLQ4kvQ74D2TnIt5Bdt4iZOcvfjN9P/B5jGZNM0flrjKqzhiX3lsi6YH0mauVop2kAyWtS7Og6yS9se/vWu7XGMqhwB2S7ic7zHVdRHwL+DhwSTp38SCycxhJXw9K9y8BXvPHMmu6qqZqcxnjlgPHAGdLOqbrsaIZyk7GuI/1KHoV2Zjh4nQtS/c/AdyeZkFvp8S/v9pWmEbE/WRTvd33HyXLhtV9v/A8RrNxIGBedd2W3RnjACR1Msbl000WzVC+BPxA0lun1S87If2AiPhher0GOAv4dirr3enRLwN3kv2HvpCXp5tVaIDY0S/pU6+McSd2lVFqhrLr+c1dZXZmOw+JiC2prC2SDu73C7QmeJiNmgZbet4v6VMdGeOGyjBXxEmfzCpU4fL0QTLGdVJI9puh3Mz0xNb5Mp/qJH5KX5/uV8HWtDy8wtSaoMJ1HrszxgFPkiV0+mDXM50Zyh9SYoYydUdeTNtE7gbOBf53V1mXM30WtFBrgofZqInqDgOqI2NcRPwY+DOyTa+vIxso/Xb6yOXATZI+DDxBickLBw+zqgywhqOMOjLGRcR64J097j8DnDpI/Rw8zCqkCTrFtDXBw2MeNmpOvWBmQ3PwMLOhTNLGuNYED++qtVHLUi+MuhazpzXBw6wJ2nK4cRkOHmYV8YCpmQ1tghoe7Qke3lVroyfmeJ2HmQ1KuOVhZsMQzJugQY/WBI9X3W2xEXPLw8yG5qlaMxvKBMUOBw+zqojJOpqvNcFj165do66CTTpNVrdlkgKlWa2yFaYqdZUqb8ikT+m9S9P9n0j6j+ne2yXdl7tekHRxeu8ySU/m3ju9X/1a0/Iwa4Kq2h25pE+nkR1cfI+ktekowY7dSZ8krSBL+vSBlBxqBfCbwK8D35P0toj4CXBcrvwnyXJId1wZEZ8vW8fWBA9P1VoTVNhrGTrpU7p/Q0RsB36Wzjh9F9lByR2nAj+NiMeHraC7LWaVEVK5i5T0KXet7CqsV9KnNxU9kxLHd5I+lfnsCuD6rnsXSrpf0nUjzVUr6fCUbHeTpAclXZTu90yoq8zVqZ92v6QT6qqbWR06sy1lLlLSp9y1ukdx3comfZrxs5L2Ac4A/k/u/VXAW8i6NVuAv+xRxjR1tjymgI9GxDuAk4ALUl+sKKHucvYk311J9suYjZUKB0z3JulTv88uBzZGxFOdGxHxVETsiohXgS/SI590tzoTXW8hi2BExIuSNpE1nYoS6p4JrElJa+6StEDSoZ38mf14zMNGTpUeQzh00idJa4GvSbqCbMB0MfDPuc+dTVeXpevf2u8BP+pXwVkZME1TSMeTZakqSqhb1E+bFjxS33AlwOFHHFFrvc0GUeUisb1J+pSeu4lscHUKuCAidgFI+jWyGZw/6fqRfyHpOLLuzWM93n+N2oOHpNcDXwcujogXZojMpZLwpr7haoAlS5b64FJrlCoPQN7LpE+fBT7b4/7LZIOq3ffPGbR+tQYPSfPJAsdXI+Ib6fZTnSZSV0LdMn28Ql5hak0wOetL651tEVmzalNEXJF7q9NPg+kJddcC56ZZl5OA58uOd5g1gYC5UqmrDepseZwMnAM8IOm+dO+TFCfUvRU4HXgEeBk4v8a6mdWiJXGhlDpnW35AcSvuNQl10yzLBXXVx6x+cq7aceSpWmsCtzzMbGDZVO3kRA8HD7OqyC2PseRuizXBJB0G1JrgYTZqTjdpZkPzbIuZDWWCei3tCR5enm5N4JaHmQ3MYx5mNpwBTkZvg9YED0/VWhNMTuhoUfAwG7VO3pZJ4dPTzSqkklepsipO+pTuPybpgZTYaX3ufs+DyWfi4GFWpYqiRy7p03LgGODsdIB43u6kT8CVZEmf6Er6tAz4Qiqv47cj4riIWJq7V3QweaHWdFs85mFNUGG3pe6kT92KDiYv5JaHWYUGaHiMKulTAN+VtKHrZ047mBw4mD5a0/Iwa4TyDY9tXd2GMiVVkfTp5Ij4RcpasE7SQxHx/f7Vfa3WBI9XX3W3xUYra1VU1m0ZJOnT5rJJnyKi8/VpSbeQdWe+T/HB5IXcbTGrSjrPo8xVwu6kTyk95AqyQ8Lz8oeJ7076lO6vSLMxR5OSPknaX9IbACTtD7yHPcmdig4mL9SalodZE1TV7qgj6ZOkQ4BbUm6ZecDXIuI76UcWHUxeyMHDrDJqdNKnNHNzbMHzz9DjYPKZtCd47JoadQ3MvCXfzAY3yOrRNnDwMKvSBEWP9gSPXTtHXQMzHwZkZsOZpDGPOhNdXyfpaUk/yt3ruXMvJbe+Ou0CvF/SCXXVy6w21a7zaLw6F4n9LdmOvryinXvLyRayLAZWAqtqrJdZbVTyf21QZ6Lr7+fPF0iKdu6dCaxJq+PukrSgs1S29A/0mIeNmGhPq6KM2V6eXrRzr8wOQgAkrezsRNy6bWutlTUbVJWHATVdU/a2lNlBmN2MWB0RSyNi6aKFi2qultmAJih6zPZsS9HOvTI7CGf2qvO22Oj5DNP6FO3cWwucm2ZdTgKeH2i8w6whJqjhUV/LQ9L1ZIOjCyVtBv47xTv3bgVOBx4BXgbOr6teZrVqS2Qooc7ZlrML3nrNzr00y3JBXXUxmw0VHwbUeO1ZYepdtTZqLVoAVkZ7godZA0xQ7GjMVK1ZC2SHAZW5SpVWcdInSYdLukPSJkkPSroo9/xlkp5MyaDuk3R6v/q1p+XhFabWAFV1W3JJn04jW8pwj6S1EZHP27I76ZOkFWRJnz7QlfTp14HvSXob2ZGEH42Ijeks0w2S1uXKvDIiPl+2jm55mFWk7DRtyfiyO+lTROwAOkmf8s4k2+YBWdKnU7uTPkXEz8hmMd8VEVsiYiNARLwIbKJgJXcZDh5mVaouetSV9CmrZtbFOR64O3f7wrSr/TrnqjWbZQPsqu2XMa6upE9Iej3wdeDiiHgh3V4FvAU4DtgC/GW/37VFYx5enm6jN8CYR7+McbUkfZI0nyxwfDUivtF5ICKe2vM76IvAt/r9Am55mFVFMKfkVUIdSZ9ElutlU0RcMa3q2V6zjt9jTzKoQu1peZg1QjXTLTUlffot4BzgAUn3pR/1yZQf5i8kHUfWvXkM+JN+dWxP8PBUrY1Y1YcB1ZD06QcURLeIOGfQ+rUneJg1wCStMHXwMKuQ97aY2VCqzFXbdO0JHt5Vaw0wOaGjTcHDbMTalJOlDAcPswr5MKBx5AOQrQkmJ3a0KHiYNcAExQ4HD7PqaKJSLzh4mFVk0tJNtid4eHm62axqT/AwawC3PMxsKJ6qHUeeqrVR8yIxMxvGpA2YNuoksX55KsyaboAzTMdeY4JHLk/FcuAY4OyUf8JsbHT2t/S72qBJ3ZbdeSoAJHXyVPx4xk8lz95zTY1VMyunyrggaRnwV2THEF4bEZd3vb8vsAZYAjwDfCAiHkvvXUqWFGoX8JGIuG2mMtNZpzcABwIbgXNSvphCjWl5UC5PBZJWdo6r37pt66xVzqyUivK2lGyJ784YB1xJljGOroxxy4AvSJrbp8zPkWWMWww8m8qeUZOCR5k8FUTE6ohYGhFLFy1cNAvVMitHwByp1FVC5RnjispMnzkllUEq86x+FWxSt6VMnoppNm7csO118/U4sBDYVmPdxpX/Lr31+rscubeFbty44bbXzdfCko/vJ2l97vXqiFide92rJX5iVxnTMsZJymeMu6vrs51WfK8yDwKeS1nnup8v1KTgsTtPBfAkWbPrgzN9ICIWAUha3yeBzkTy36W3uv4uEbGswuLqyBjXq6fRN8NckcYEj6I8FSOultmo1JIxruD+NmCBpHmp9dG31Q/NGvMgIm6NiLdFxFtS3gmzSVV5xriiMtNn7khlkMr8Zr8KNqblsZdW939kIvnv0lvj/y51ZIwDmKF1/3HgBkmfAe5NZc9IWdAxMxtMo7otZjY+HDzMbChjHTwmeSOdpMMl3SFpk6QHJV2U7h8oaZ2kh9PXN6b7knR1+lvdL+mE0f4G9UorKu+V9K30+mhJd6e/y41pwJA0qHhj+rvcLemoUdZ7nIxt8PBGOqaAj0bEO4CTgAvS7/8J4Pa0zPj29Bqyv9PidK0EVs1+lWfVRcCm3Oui5dc9l3hbf2MbPCi3fLe1ImJLRGxM379I9g/lTUxfspxfZnwmsCYyd5HN6x86y9WeFZIOA34HuDa9nmn5ddESb+tjnINHqY10kyA1tY8H7gYOiYgtkAUY4OD02CT9va4C/hx4Nb2eafn1tCXeQGeJt/UxzsFjqCW1bSPp9cDXgYsj4oWZHu1xr3V/L0nvBZ6OiA352z0ejRLv2QzGeZHYwBvp2kbSfLLA8dWI+Ea6/ZSkQyNiS+qWPJ3uT8rf62TgDEmnA/sBB5C1RIqWXxct8bY+xrnlUWb5bmulfvmXgE0RcUXurfyS5fwy47XAuWnW5STg+U73pk0i4tKIOCwijiL7/8Q/RMSHKF5+XbTE2/oY6xWm6b8uV7Fnqe3E7IeR9FvAPwIPsKdv/0mycY+bgCOAJ4D3R8QvU7C5huxwmJeB8yNi/WsKbhFJ7wY+FhHvlfRm9pyUdS/whxGxXdJ+wFfIxox+CazonGZnMxvr4GFmozPO3RYzGyEHDzMbioOHmQ3FwcPMhuLgYWZDcfBoqbTr9meSDkyv35he7/Up4Wbg4NFaEfFzsp2znSxjl5Md7//46GplbeJ1Hi2Wlq9vAK4D/hg4vl8KQbOyxnlvi/URETsl/VfgO8B7HDisSu62tN9yYAvwzlFXxNrFwaPFJB0HnEZ20th/aevhPzYaDh4tlTbCrSI75+MJ4H8Bnx9traxNHDza64+BJyJiXXr9BeA3JP37EdbJWsSzLWY2FLc8zGwoDh5mNhQHDzMbioOHmQ3FwcPMhuLgYWZDcfAws6H8fxWGP2kgrtIuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imshow_grid(rg, 'surface_water__depth', cmap='Blues')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks pretty good. Let's check the values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.        ,  0.02153325,  0.02153325,  0.02153325,  0.02153325,\n",
       "        0.02153325,  0.02153325,  0.02153325,  0.02153325,  0.02153325,\n",
       "        0.02153325,  0.02153325,  0.02153325,  0.02153325,  0.02153325,\n",
       "        0.02153325,  0.02153325,  0.02153325,  0.02153325,  0.02153325,\n",
       "        0.02153325,  0.02153325,  0.02153325,  0.02153325,  0.02153325,\n",
       "        0.02153325,  0.02153325,  0.02153325,  0.02153325,  0.02153325,\n",
       "        0.02153325,  0.02153325,  0.02153325,  0.02153325,  0.02153325,\n",
       "        0.02153325,  0.02153325,  0.02153325,  0.02153325,  0.02153325,\n",
       "        0.02153325,  0.        ])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rg.at_node['surface_water__depth'][42:84] # bottom row of core nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the depth agrees with the analytical solution to within three decimal places: not bad. Ideally, we would build the above tests into the component as doctests or unit tests. We could also test the transient solutions: rising hydrograph, falling hydrograph. Finally, we haven't tested all the ingredients; for example, we haven't tested what happens when infiltration rate is greater than zero.\n",
    "\n",
    "Nonetheless, the above example illustrates the basics of component-making. A great next step would be to create a unit test based on this example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next: <a href=\"https://github.com/landlab/landlab/wiki/Tutorials\">Tips on Writing Landlab Components (URL to be added)</a>\n",
    "\n",
    "### Click here for more <a href=\"https://github.com/landlab/landlab/wiki/Tutorials\">Landlab tutorials</a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}