{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Profiling and Scaling Analysis of the NetworkSedimentTransporter\n", "\n", "Planning to run the NetworkSedimentTransporter (NST) on a large river network, or for a long time? You might wonder: how efficient is this component? \n", "\n", "Most field-scale applications of the NST would likley have 100-500 links in the network and be initialized with at least 100 parcels per link (e.g., up to 50k+ parcels). This number of parcels per link are commonly need to capture the full grain size distribution. \n", "\n", "Because the component may use many parcels on many links, we will want to assess two parts of the components efficiency. First, we will *profile* the code to see what happens each time the component runs. This is useful to see which parts of the computation take the most time and occur the most often. If code is running too slowly, it is usefull to profile it in order to identify which parts may yield the biggest speedups. \n", "\n", "The second goal is *scaling*, or seeing how changes in the number of parcels and links changes run time. If doubling the number of links doubled the run time, we would say that runtime for the component \"scales linearly\" with the number of links. In computer science this is often called the [time complexity](https://en.wikipedia.org/wiki/Time_complexity) of the code/algorithm. \n", "\n", "The goal of this notebook is to explore profiling and scaling with the NST. We will look at how changing the grid nodes and number of parcels per link impacts the initialization and run time of this component. In order to do this we will need to make some python functions to create variably sized grids. We will also use some common python tools for profiling code. The notebook is organized into three parts:\n", "\n", "- Part 1: Create generic, variable sized grids and parcels. \n", "\n", "- Part 2: Profiling the code.\n", "\n", "- Part 3: Scaling analysis.\n", "\n", "We begin by importaing all needed python modules. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import cProfile\n", "import io\n", "import pstats\n", "import time\n", "import warnings\n", "from pstats import SortKey\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "\n", "from landlab.components import FlowDirectorSteepest, NetworkSedimentTransporter\n", "from landlab.data_record import DataRecord\n", "from landlab.grid.network import NetworkModelGrid\n", "from landlab.plot import graph\n", "\n", "warnings.filterwarnings('ignore')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 1: Variably-sized grids and parcels\n", "\n", "### Part 1a: Grid topology\n", "\n", "First, we need the ability to create different sizes of network model grids. We will use these different size grids to assess how long it takes the NST to run with different grid sizes. \n", "\n", "A simple approach is to create a generic grid in which each node has two recievers. Lets start by writing a function that creates the x and y node coordinates and the linking structure for a given number of layers. \n", "\n", "While these grids do not look like typical river networks, they are topologically similar. \n", "\n", "We create a function `create_node_xy_and_links` which takes the following inputs\n", "\n", "- `n_layers`\n", "- `x0`\n", "- `y0`\n", "- `x_perc`\n", "- `dy`\n", "\n", "and returns three arrays:\n", "- `x_of_node`\n", "- `y_of_node`\n", "- `nodes_at_link`\n", "\n", "These inputs and outputs are defined in the function docstring below. The function was designed to produce output that could be directly provided to the [`NetworkModelGrid`](https://landlab.readthedocs.io/en/master/reference/grid/network_api.html#module-landlab.grid.network) init function. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def create_node_xy_and_links(n_layers, x0=0., y0=0., xperc=0.9, dy=1.):\n", " \"\"\"Create node and link structure of a branching binary tree.\n", " \n", " The tree can have an arbitrary number of layers. For example, \n", " a tree with one layer has three nodes and two links:\n", " \n", " ::\n", " \n", " * *\n", " \\ /\n", " *\n", " \n", " The lowest of the nodes is the \"origin\" node, and has coordinates\n", " of `(x0, y0)`. The y spacing between layers is given by `dy`. Finally, \n", " in order to ensure that links do not cross and nodes are not co-located\n", " a shrinking factor, `xperc` that must be less than 1.0 is specified. \n", " \n", " Each layer has 2^{layer} nodes: Layer 0 has 1 node, layer 1 has 2 nodes, \n", " layer 2 has 4 nodes. \n", " \n", " A tree with three layers has seven nodes and six links:\n", " \n", " ::\n", " \n", " * * * *\n", " \\ / \\ /\n", " * *\n", " \\ /\n", " *\n", " \n", " Parameters\n", " ----------\n", " n_layers : int\n", " Number of layers of the binary tree\n", " x0 : float\n", " x coordinate position of the origin node. Default of 0.\n", " y0=0. : float\n", " y coordinate position of the origin node. Default of 0. \n", " xperc : float\n", " x direction shrinkage factor to prevent co-location of nodes\n", " and corssing links. Must be between 0.0 and 1.0 noninclusive. \n", " Default of 0.9.\n", " dy : float\n", " y direction spacing between layers. Default of 1.\n", " \n", " Returns\n", " ------\n", " x_of_node : list\n", " Node x coordinates.\n", " y_of_node : list\n", " Node y coordinates.\n", " nodes_at_link : list of (head, tail) tuples\n", " Nodes at link tail and head.\n", " \n", " \"\"\"\n", " assert xperc<1.0\n", " assert xperc>0.0\n", " nodes_per_layer = np.power(2, np.arange(n_layers+1))\n", " nnodes = np.sum(nodes_per_layer) \n", " x_of_node=[x0]\n", " y_of_node=[y0]\n", " nodes_at_link = []\n", " id_start_layer = 0\n", " for nl in np.arange(1, n_layers+1):\n", " nodes_last_layer = np.power(2, nl-1)\n", " nodes_this_layer = np.power(2, nl) \n", " dx = xperc * (dy)*(0.5**(nl-1))\n", " for ni in range(nodes_last_layer):\n", " head_id = id_start_layer+ni \n", " tail_id = len(x_of_node) \n", " x = x_of_node[head_id]\n", " y = y_of_node[head_id] \n", " x_of_node.extend([x-dx, x+dx])\n", " y_of_node.extend([y+dy, y+dy]) \n", " nodes_at_link.extend([(head_id, tail_id), (head_id, tail_id +1)]) \n", " id_start_layer = len(x_of_node) - nodes_this_layer\n", " return x_of_node, y_of_node, nodes_at_link" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's demonstrate the different sorts of grids we get with different numbers of layers. We'll look at grids with between 3 and 1023 nodes. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "example_layers = [1, 3, 5, 7, 9]\n", "\n", "nodes = []\n", "for i, n_layers in enumerate(example_layers):\n", " x_of_node, y_of_node, nodes_at_link = create_node_xy_and_links(n_layers)\n", " grid = NetworkModelGrid((y_of_node, x_of_node), nodes_at_link)\n", "\n", " graph.plot_graph(grid, at=\"node,link\", with_id=False)\n", " nodes.append(grid.number_of_nodes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see the number of nodes grows exponentially with the number of layers. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(example_layers, nodes)\n", "plt.xlabel(\"Number of Layers\")\n", "plt.ylabel(\"Number of Nodes\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1b: Generic grid.\n", "\n", "In order to run the NST, our grids need fields related to channel geometry and flow characteristics. Here, we'll create a function that takes only the number of layers and uses the function `create_node_xy_and_links` that we just created, creates a grid, adds, those fields and populate them with generic values. The NST also needs a [`FlowDirectorSteepest`](https://landlab.readthedocs.io/en/master/reference/components/flow_director.html?highlight=FlowDirectorSteepest#module-landlab.components.flow_director.flow_director_steepest) instance, so we'll create that too. \n", "\n", "Thus, with these two functions we can specifiy the desired number of layers and get both of these objects we need to instantiate the NST. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def create_nmg_and_fd(n_layers):\n", " \"\"\"Create a generic NetworkModelGrid and FlowDirectorSteepest. \n", " \n", " This function will also add the following fields to the NetworkModelGrid.\n", " \n", " - topographic__elevation at node\n", " - bedrock__elevation at node\n", " - flow_depth at link\n", " - reach_length at link\n", " - channel_width at link. \n", " \n", " Parameters\n", " ----------\n", " n_layers : int\n", " Number of layers of the binary tree\n", " \n", " Returns\n", " -------\n", " grid : NetworkModelGrid instance\n", " fd : FlowDirectorSteepest instance\n", " \"\"\"\n", " x_of_node, y_of_node, nodes_at_link = create_node_xy_and_links(n_layers)\n", " grid = NetworkModelGrid((y_of_node, x_of_node), nodes_at_link)\n", "\n", " _ = grid.add_field(\"topographic__elevation\", grid.y_of_node.copy(), at=\"node\")\n", " _ = grid.add_field(\"bedrock__elevation\", grid.y_of_node.copy(), at=\"node\")\n", " _ = grid.add_field(\"flow_depth\", 2.5*np.ones(grid.number_of_links), at=\"link\") # m\n", " _ = grid.add_field(\"reach_length\", 200.*np.ones(grid.number_of_links), at=\"link\") # m\n", " _ = grid.add_field(\"channel_width\", 1.*np.ones(grid.number_of_links), at=\"link\") # m\n", " \n", " fd = FlowDirectorSteepest(grid)\n", " fd.run_one_step()\n", " return grid, fd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1c: Generic sets of parcels\n", "Our second to last step in writing functions to create a generic model grid and set of parcels is to create a function that will create our parcels. The parcels are associated with the grid because parcels must be located on a grid element, and thus we will always need the grid to create the parcels. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def create_parcels(grid, parcels_per_link=5):\n", " \"\"\"Create a generic set of parcels.\n", " \n", " The NetworkSedimentTransporter requires a set of parcels with some \n", " specific attributes (e.g., density) that are used in order to \n", " calculate travel distances. This function creates the parcels in \n", " the correct format and populates all necessary attributes. Specifically\n", " it creates the following attributes:\n", " \n", " - \"abrasion_rate\"\n", " - \"density\"\n", " - \"time_arrival_in_link\"\n", " - \"active_layer\"\n", " - \"location_in_link\"\n", " - \"D\"\n", " - \"volume\"\n", " \n", " Parameters\n", " ----------\n", " grid : NetworkModelGrid\n", " parcels_per_link : int\n", " Number of parcels to create for each link. Default of 5. \n", " \n", " Returns\n", " -------\n", " parcels : DataRecord\n", " \n", " \"\"\"\n", " # element_id is the link on which the parcel begins. \n", " element_id = np.repeat(np.arange(grid.number_of_links), parcels_per_link)\n", " element_id = np.expand_dims(element_id, axis=1)\n", "\n", " # scale volume with parcels per link so we end up with a similar quantity of sediment. \n", " volume = (1./parcels_per_link) * np.ones(np.shape(element_id)) # (m3)\n", " \n", " active_layer = np.zeros(np.shape(element_id)) # 1= active, 0 = inactive\n", " density = 2650 * np.ones(np.size(element_id)) # (kg/m3)\n", " abrasion_rate = 0. * np.ones(np.size(element_id)) # (mass loss /m)\n", "\n", " # Lognormal GSD\n", " medianD = 0.085 # m\n", " mu = np.log(medianD)\n", " sigma = np.log(2) #assume that D84 = sigma*D50\n", " np.random.seed(0)\n", " D = np.random.lognormal(\n", " mu,\n", " sigma,\n", " np.shape(element_id)\n", " ) # (m) the diameter of grains in each parcel\n", "\n", " time_arrival_in_link = np.random.rand(np.size(element_id), 1) \n", " location_in_link = np.random.rand(np.size(element_id), 1) \n", "\n", " variables = {\n", " \"abrasion_rate\": ([\"item_id\"], abrasion_rate),\n", " \"density\": ([\"item_id\"], density),\n", " \"time_arrival_in_link\": ([\"item_id\", \"time\"], time_arrival_in_link),\n", " \"active_layer\": ([\"item_id\", \"time\"], active_layer),\n", " \"location_in_link\": ([\"item_id\", \"time\"], location_in_link),\n", " \"D\": ([\"item_id\", \"time\"], D),\n", " \"volume\": ([\"item_id\", \"time\"], volume)\n", " }\n", "\n", "\n", " items = {\"grid_element\": \"link\", \"element_id\": element_id}\n", "\n", " parcels = DataRecord(\n", " grid,\n", " items=items,\n", " time=[0.0],\n", " data_vars=variables,\n", " dummy_elements={\"link\": [NetworkSedimentTransporter.OUT_OF_NETWORK]},\n", " )\n", "\n", " return parcels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2: Profile the code\n", "\n", "Next we will create a NST, run it, and profile it using the python built-in profiling tool [cProfile](https://docs.python.org/3.8/library/profile.html).\n", "\n", "We create our grid and parcels, instantiate the profile, run the model for 60 timesteps, and then compile the results of profiling. \n", "\n", "We will profile the code in two parts, the instantiation and running. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# feel free to change these parameters and see\n", "# how it impacts the results\n", "nlayer = 5\n", "timesteps = 50\n", "parcels_per_link = 50\n", "\n", "# calculate dt and set seed. \n", "dt = 60 * 60 * 24 *12 # length of timestep (seconds) \n", "np.random.seed(1234)\n", "\n", "pr = cProfile.Profile()\n", "pr.enable()\n", "\n", "grid, fd = create_nmg_and_fd(nlayer)\n", "\n", "parcels = create_parcels(grid, parcels_per_link=parcels_per_link)\n", "\n", "nst = NetworkSedimentTransporter( \n", " grid,\n", " parcels,\n", " fd,\n", " bed_porosity=0.3,\n", " g=9.81,\n", " fluid_density=1000,\n", " transport_method=\"WilcockCrowe\",\n", " )\n", "\n", "pr.disable()\n", "s = io.StringIO()\n", "sortby = SortKey.CUMULATIVE\n", "ps = pstats.Stats(pr, stream=s).sort_stats(sortby)\n", "ps.print_stats()\n", "print(s.getvalue())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the majority of the time it took to create the grid, parcels, and the NST component occurs in `xarray` functions. This is good, as those are likely already efficient. \n", "\n", "Next we profile the code as we use the `run_one_step` function that advances the component forward in time. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pr = cProfile.Profile()\n", "pr.enable()\n", "for t in range(0, (timesteps * dt), dt):\n", " nst.run_one_step(dt)\n", "pr.disable()\n", "s = io.StringIO()\n", "sortby = SortKey.CUMULATIVE\n", "ps = pstats.Stats(pr, stream=s).sort_stats(sortby)\n", "ps.print_stats()\n", "print(s.getvalue())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the majority of the time it takes to run the component happens in a function called `_partition_active_and_storage_layers`. This function is used to figure out which parcels are moving and which are not active. \n", "\n", "\n", "## Part 3: Scaling\n", "Next we look at how initialization and moving the model forward in time change as we increase the number of nodes and number of parcels on a given link. \n", "\n", "### 3a: Create a timing function. \n", "To do this we create one final convenience function that will create everything we need for our scaling analysis, run the model, and report how long it took to initialize and run. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def time_code(nlayer=3, parcels_per_link=100, timesteps=10):\n", " \"\"\"Time the initializiation and runtime.\n", " \n", " Parameters\n", " ----------\n", " n_layers : int\n", " Number of layers of the binary tree used to create the \n", " NetworkModelGrid. Default of 3\n", " parcels_per_link : int\n", " Number of parcels to create for each link. Default of 100.\n", " timesteps : int\n", " Number of timesteps. Default of 10. \n", " \n", " Returns \n", " -------\n", " (number_of_nodes, parcels_per_link, total_parcels) : tuple\n", " Tuple indicating the key inputs in our scaling analysis. The number \n", " of nodes, the number of parcels per link, the total number of parcels. \n", " init_duration : float\n", " Duration of the initiation step, in seconds.\n", " r1s_per : float\n", " Duration of the average run_one_step call, in seconds\n", " \"\"\"\n", " init_start = time.time()\n", "\n", " grid, fd = create_nmg_and_fd(nlayer)\n", "\n", " parcels = create_parcels(grid, parcels_per_link=parcels_per_link)\n", " \n", " dt = 60 * 60 * 24 *12 # length of timestep (seconds) \n", "\n", " nst = NetworkSedimentTransporter( \n", " grid,\n", " parcels,\n", " fd,\n", " bed_porosity=0.3,\n", " g=9.81,\n", " fluid_density=1000,\n", " transport_method=\"WilcockCrowe\",\n", " )\n", " init_duration = time.time() - init_start\n", "\n", " if timesteps >0:\n", " r1s_start = time.time()\n", " for t in range(timesteps):\n", " nst.run_one_step(dt)\n", " r1s_per = (time.time() - r1s_start) / timesteps\n", " else: \n", " r1s_per = 0.0\n", "\n", " return (grid.number_of_nodes, parcels_per_link), init_duration, r1s_per" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we use or new `time_code` function with a few different grid sizes, a few different parcels per link, and for 10 seconds. Feel free to experiment and change these values. Some of these values have been reduced to ensure that this notebook always works in the Landlab continuous integration. \n", "\n", "### 3b: Time the code while systematically varying inputs\n", "We save the results in a list structured to easily make a [`pandas.Dataframe`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.random.seed(345)\n", "\n", "out = []\n", "# this range for i in reduced for testing. \n", "for i in range(2, 5):\n", " for j in [10, 20, 50, 100, 200, 500]:\n", " print(i, j)\n", " (nn, ppl), init, r1s_per = time_code(\n", " nlayer=i, \n", " parcels_per_link=j, \n", " timesteps=10)\n", " out.append((nn, ppl, init, r1s_per))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We make a dataframe and investigate the contents with `df.head`. We'll use some shorthand for the column and axis names:\n", "\n", "- nnodes = number of nodes in the NetworkModelGrid\n", "- ppl = parcels per link\n", "- init = duration of time (in seconds) of instantiating the grid, parcels, and NST\n", "- r1s_per = duration of time (in seconds) of running the model forward one step in time. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "df = pd.DataFrame(out, columns=[\"nnodes\", \"ppl\", \"init\", 'r1s_per'])\n", "df = df.pivot(index='nnodes', columns='ppl', values=[\"init\", 'r1s_per'])\n", "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3c: Plot the output to see how the code scales. \n", "\n", "First, lets plot the duration of init (left) and run one step (right) as a function of the number of nodes and parcels per link. We will put number of nodes on the x axis and make a line for each different nubmer of parcels per link. For each As is common, we will make these plots in log-log space." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, dpi=300)\n", "\n", "df[\"init\"].plot(loglog=True, ax=ax[0], title=\"init duration\")\n", "ax[0].set_ylabel(\"duration (s)\")\n", "ax[0].set_xlabel(\"Number of Nodes\")\n", "\n", "df[\"r1s_per\"].plot(loglog=True, ax=ax[1], title=\"run one step duration\")\n", "ax[1].set_ylabel(\"duration (s)\")\n", "ax[1].set_xlabel(\"Number of Nodes\")\n", "#plt.savefig(\"scaling1.png\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that both the init and run one step functions increase in their duration when the number of nodes parcels per link increase. However, the init function is less sensitive to the number of parcels per link. Feel free to explore plotting this data in different ways. \n", "\n", "### Click here for more Landlab tutorials" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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": 4 }