{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Mapping values between grid elements\n", "\n", "
\n", "For more Landlab tutorials, click here: https://landlab.readthedocs.io/en/latest/user_guide/tutorials.html)\n", "
\n", "\n", "Imagine that you're using Landlab to write a model of shallow water flow over terrain. A natural approach is to place your scalar fields, such as water depth, at the nodes. You then place your vector fields, such as water surface gradient, flow velocity, and discharge, at the links. But your velocity depends on both slope and depth, which means you need to know the depth at the links too. How do you do this?\n", "\n", "This tutorial introduces *mappers*: grid functions that map quantities defined on one set of elements (such as nodes) onto another set of elements (such as links). As you'll see, there are a variety of mappers available.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mapping from nodes to links\n", "\n", "For the sake of example, we'll start with a simple 3-row by 4-column raster grid. The grid will contain a scalar field called `water__depth`, abbreviated `h`. We'll populate it with some example values, as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from landlab import RasterModelGrid\n", "import numpy as np\n", "mg = RasterModelGrid((3, 4), xy_spacing=100.0)\n", "h = mg.add_zeros('surface_water__depth', at='node')\n", "h[:] = 7 - np.abs(6 - np.arange(12))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the sake of visualizing values at nodes on our grid, we'll define a handy little function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def show_node_values(mg, u):\n", " for r in range(mg.number_of_node_rows - 1, -1, -1):\n", " for c in range(mg.number_of_node_columns):\n", " print(int(u[c + (mg.number_of_node_columns * r)]), end=' ')\n", " print()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "show_node_values(mg, h)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's review the numbering of nodes and links. The lines below will print a list that shows, for each link ID, the IDs of the nodes at the link's tail and head:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i in range(mg.number_of_links):\n", " print(i, mg.node_at_link_tail[i], mg.node_at_link_head[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Finding the mean value between two nodes on a link\n", "\n", "Suppose we want to have a *link-based* array, called *h_edge*, that contains water depth at locations between adjacent pairs of nodes. For each link, we'll simply take the average of the depth at the link's two nodes. To accomplish this, we can use the `map_mean_of_link_nodes_to_link` grid method. At link 8, for example, we'll average the *h* values at nodes 5 and 6, which should give us a depth of (6 + 7) / 2 = 6.5:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_edge = mg.map_mean_of_link_nodes_to_link('surface_water__depth')\n", "for i in range(mg.number_of_links):\n", " print(i, h_edge[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What's in a name?\n", "\n", "The mapping functions have long names, which are designed to make it as clear as possible to understand what each function does. All the mappers start with the verb *map*. Then the relationship is given; in this case, we are looking at the *mean*. Then the elements from which a quantity is being mapped: we are taking values from *link nodes*. Finally, the element to which the new values apply: *link*.\n", "\n", "### Mapping minimum or maximum values\n", "\n", "We can also map the minimum value of *h*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_edge = mg.map_min_of_link_nodes_to_link('surface_water__depth')\n", "for i in range(mg.number_of_links):\n", " print(i, h_edge[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... or the maximum:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_edge = mg.map_max_of_link_nodes_to_link('surface_water__depth')\n", "for i in range(mg.number_of_links):\n", " print(i, h_edge[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Upwind and downwind\n", "\n", "Numerical schemes often use *upwind differencing* or *downwind differencing*. For example, finite difference schemes for equations that include advection may use \"upwind\" rather than centered differences, in which a scalar quantity (our *h* for example) is taken from whichever side is upstream in the flow field.\n", "\n", "How do we know the flow direction? If the flow is driven by the gradient in some scalar field, such as pressure or elevation, one approach is to look at the values of this scalar on either end of the link: the end with the higher value is upwind, and the end with the lower value is downwind.\n", "\n", "Suppose for example that our water flow is driven by the water-surface slope (which is often a good approximation for the *energy slope*, though it omits the kinetic energy). Let's define a bed-surface elevation field *z*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "z = mg.add_zeros('topographic__elevation', at='node')\n", "z[:] = 16 - np.abs(7 - np.arange(12))\n", "show_node_values(mg, z)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The water-surface elevation is then the sum of *h* and *z*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "w = z + h\n", "show_node_values(mg, w)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For every link, we can assign the value of *h* from whichever end of the link has the greater *w*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_edge = mg.map_value_at_max_node_to_link(w, h)\n", "for i in range(mg.number_of_links):\n", " print(i, h_edge[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the middle two nodes (5 and 6). Node 6 is higher (22 versus 20). Therefore, the link between them (link 8) should be assigned the value of *h* at node 6. This value happens to be 7.0.\n", "\n", "Of course, we could also take the value from the *lower* of the two nodes, which gives link 8 a value of 6.0:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_edge = mg.map_value_at_min_node_to_link(w, h)\n", "for i in range(mg.number_of_links):\n", " print(i, h_edge[i])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### Heads or tails?\n", "\n", "It is also possible to map the scalar quantity at either the head node or the tail node to the link:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_edge = mg.map_link_head_node_to_link('surface_water__depth')\n", "for i in range(mg.number_of_links):\n", " print(i, h_edge[i])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_edge = mg.map_link_tail_node_to_link('surface_water__depth')\n", "for i in range(mg.number_of_links):\n", " print(i, h_edge[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simple example using centered water depth\n", "The following implements one time-step of a linear-viscous flow model, in which flow velocity is calculated at the links, and the depth at each link is taken as the mean of depth at the two bounding nodes. To make the flow a little tamer, we'll have our fluid be hot, low viscosity basaltic lava instead of water, with a dynamic viscosity of 100 Pa s." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gamma = 25000.0 # unit weight of fluid, N/m2\n", "viscosity = 100.0 # dynamic viscosity in Pa s\n", "grad = mg.calc_grad_at_link(w)\n", "h_edge = mg.map_mean_of_link_nodes_to_link(h)\n", "vel = -(gamma / (3.0 * viscosity)) * h_edge * h_edge * grad\n", "for ln in range(mg.number_of_links):\n", " print(ln, h_edge[ln], grad[ln], vel[ln])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I'm not sure I love the idea of a 5-m thick lava flow moving at over 100 m/s! (I guess we can take some comfort from the thought that turbulence would probably slow it down)\n", "\n", "How different would the numerical solution be using an upwind scheme for flow depth? Let's find out:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_edge = mg.map_value_at_max_node_to_link(w, h)\n", "vel = -(gamma / (3.0 * viscosity)) * h_edge * h_edge * grad\n", "for ln in range(mg.number_of_links):\n", " print(ln, h_edge[ln], grad[ln], vel[ln])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even scarier.\n", "\n", "In any event, this example illustrates how you can use Landlab's mapping functions to build mass-conservation models in which the flow rate depends on a gradient and a scalar, both of which can be evaluated at links." ] } ], "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 }