From dfea69c2769ee1dd164ae306c4d66142f6a51b8c Mon Sep 17 00:00:00 2001 From: angelmons Date: Sun, 10 Nov 2024 23:54:20 -0300 Subject: [PATCH] Update tests and rename using CamelCase Added several tests, including analytical solutions, renamed the component to RiverFlowDynamics, changes the documentation in main component file, and created a second tutorial --- .../river_flow_dynamics_tutorial.ipynb | 14 +- .../river_flow_dynamics_tutorial2.ipynb | 425 ++++++++++++++++++ notebooks/teaching | 2 +- notebooks/tutorials | 2 +- src/landlab/components/__init__.py | 4 +- .../river_flow_dynamics/__init__.py | 4 +- .../river_flow_dynamics.py | 266 ++--------- .../test_river_flow_dynamics.py | 393 +++++++++++++++- 8 files changed, 859 insertions(+), 251 deletions(-) create mode 100644 docs/source/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial2.ipynb diff --git a/docs/source/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial.ipynb b/docs/source/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial.ipynb index c3e8f6094c..8241cf03ba 100644 --- a/docs/source/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial.ipynb +++ b/docs/source/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial.ipynb @@ -68,7 +68,7 @@ "from tqdm import trange\n", "\n", "from landlab import RasterModelGrid\n", - "from landlab.components import river_flow_dynamics\n", + "from landlab.components import RiverFlowDynamics\n", "from landlab.io import esri_ascii" ] }, @@ -87,7 +87,7 @@ "metadata": {}, "outputs": [], "source": [ - "help(river_flow_dynamics)" + "help(RiverFlowDynamics)" ] }, { @@ -116,7 +116,7 @@ "channel_slope = 0.01 # Channel slope [m/m]\n", "\n", "# Simulation parameters\n", - "n_timesteps = 100 # Number of timesteps\n", + "n_timesteps = 1000 # Number of timesteps\n", "dt = 0.1 # Timestep duration, [s]\n", "nrows = 20 # Number of node rows\n", "ncols = 60 # Number of node cols\n", @@ -264,7 +264,7 @@ "outputs": [], "source": [ "# Finally, we run the model and let the water fill our channel\n", - "rfd = river_flow_dynamics(\n", + "rfd = RiverFlowDynamics(\n", " grid,\n", " dt=dt,\n", " mannings_n=mannings_n,\n", @@ -557,7 +557,7 @@ "outputs": [], "source": [ "# Finally, we run the model and let the water fill our channel\n", - "rfd = river_flow_dynamics(\n", + "rfd = RiverFlowDynamics(\n", " grid,\n", " dt=dt,\n", " mannings_n=mannings_n,\n", @@ -618,7 +618,7 @@ "-- --\n", "### And that's it! \n", "\n", - "Nice work completing this tutorial. You know now how to use the `river_flow_dynamics` Landlab component to run your own simulations :)\n", + "Nice work completing this tutorial. You know now how to use the `RiverFlowDynamics` Landlab component to run your own simulations :)\n", "\n", "-- --\n", "\n" @@ -648,7 +648,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.0" + "version": "3.12.5" } }, "nbformat": 4, diff --git a/docs/source/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial2.ipynb b/docs/source/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial2.ipynb new file mode 100644 index 0000000000..9bca504c04 --- /dev/null +++ b/docs/source/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial2.ipynb @@ -0,0 +1,425 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 2D Surface Water Flow component\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# River Flow Dynamics Simulation with Landlab" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "For more Landlab tutorials, click here: https://landlab.readthedocs.io/en/latest/user_guide/tutorials.html\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "This notebook demonstrate the usage of the `river flow dynamics` Landlab component. The component runs a semi-implicit, semi-Lagrangian finite-volume approximation to the depth-averaged 2D shallow-water equations of Casulli and Cheng (1992) and related work.\n", + "\n", + "This notebook demonstrates how to simulate river flow dynamics using the Landlab library, implementing the semi-implicit, semi-Lagrangian finite-volume approximation of the depth-averaged shallow water equations (Casulli and Cheng, 1992).\n", + "\n", + "\n", + "## Setup and Imports\n", + "\n", + "Import the needed libraries:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from landlab import RasterModelGrid\n", + "from landlab.components import RiverFlowDynamics # Note: Using updated CamelCase naming\n", + "from landlab.io import read_esri_ascii\n", + "from landlab.plot.imshow import imshow_grid" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create Grid and Set Initial Conditions\n", + "\n", + "First, let's create a rectangular grid for our flow dynamics calculations:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nRows = 20\n", + "nCols = 60\n", + "cellSize = 0.1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Creating the grid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "grid = RasterModelGrid((nRows, nCols), xy_spacing=(cellSize, cellSize))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setting up the initial topographic elevation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "te = grid.add_zeros(\"topographic__elevation\", at=\"node\")\n", + "te += 0.059 - 0.01 * grid.x_of_node\n", + "te[grid.y_of_node > 1.5] = 1.0\n", + "te[grid.y_of_node < 0.5] = 1.0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualizing the initial topography" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(12, 4))\n", + "imshow_grid(grid, \"topographic__elevation\")\n", + "plt.title(\"Initial Topographic Elevation\")\n", + "plt.colorbar(label=\"Elevation (m)\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualizing the middle bed profile" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "middleBedProfile = np.reshape(te, (nRows, nCols))[10, :]\n", + "plt.figure(figsize=(12, 3))\n", + "plt.plot(middleBedProfile)\n", + "plt.title(\"Middle Longitudinal Section of Bed Profile\")\n", + "plt.xlabel(\"Distance (cells)\")\n", + "plt.ylabel(\"Elevation (m)\")\n", + "plt.grid(True)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Initializing Required Fields\n", + "\n", + "Create water depth field (initially empty channel)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "h = grid.add_zeros(\"surface_water__depth\", at=\"node\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create velocity field (initially zero)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "vel = grid.add_zeros(\"surface_water__velocity\", at=\"link\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Calculate initial water surface elevation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "wse = grid.add_zeros(\"surface_water__elevation\", at=\"node\")\n", + "wse += h + te" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setting up the boundary conditions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fixed_entry_nodes = np.arange(300, 910, 60)\n", + "fixed_entry_links = grid.links_at_node[fixed_entry_nodes][:, 0]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Set fixed values for entry nodes/links" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "entry_nodes_h_values = np.full(11, 0.5) # 0.5m water depth\n", + "entry_links_vel_values = np.full(11, 0.45) # 0.45 m/s velocity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run Simulation\n", + "\n", + "Initialize the RiverFlowDynamics component" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rfd = RiverFlowDynamics(\n", + " grid,\n", + " dt=0.1,\n", + " mannings_n=0.012,\n", + " fixed_entry_nodes=fixed_entry_nodes,\n", + " fixed_entry_links=fixed_entry_links,\n", + " entry_nodes_h_values=entry_nodes_h_values,\n", + " entry_links_vel_values=entry_links_vel_values,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Run the simulation for 100 timesteps (10 seconds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n_timesteps = 100\n", + "for timestep in range(n_timesteps):\n", + " rfd.run_one_step()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Analyze Results\n", + "\n", + "Get flow depth along center of channel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "flow_depth = np.reshape(grid[\"node\"][\"surface_water__depth\"], (nRows, nCols))[10, :]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot flow depth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(12, 4))\n", + "plt.plot(flow_depth, label='Simulated')\n", + "plt.title(\"Flow Depth Along Channel Centerline\")\n", + "plt.xlabel(\"Distance (cells)\")\n", + "plt.ylabel(\"Depth (m)\")\n", + "plt.grid(True)\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Get and plot velocity along center of channel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "linksAtCenter = grid.links_at_node[np.array(np.arange(600, 660))][:-1, 0]\n", + "flow_velocity = grid[\"link\"][\"surface_water__velocity\"][linksAtCenter]\n", + "\n", + "plt.figure(figsize=(12, 4))\n", + "plt.plot(flow_velocity, label='Simulated')\n", + "plt.title(\"Flow Velocity Along Channel Centerline\")\n", + "plt.xlabel(\"Distance (cells)\")\n", + "plt.ylabel(\"Velocity (m/s)\")\n", + "plt.grid(True)\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization of Final State\n", + "\n", + "Create a figure with two subplots and then let's plot final water depth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))\n", + "\n", + "# Plot final water depth\n", + "plt.subplot(2, 1, 1)\n", + "im1 = imshow_grid(grid, \"surface_water__depth\")\n", + "plt.title(\"Final Water Depth\")\n", + "plt.colorbar(label=\"Depth (m)\")\n", + "\n", + "# Plot final water surface elevation\n", + "plt.subplot(2, 1, 2)\n", + "im2 = imshow_grid(grid, \"surface_water__elevation\")\n", + "plt.title(\"Final Water Surface Elevation\")\n", + "plt.colorbar(label=\"Elevation (m)\")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "-- --\n", + "### And that's it! \n", + "\n", + "Nice work completing this tutorial. You know now how to use the `RiverFlowDynamics` Landlab component to run your own simulations :)\n", + "\n", + "-- --\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Click here for more Landlab tutorials" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.12.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/teaching b/notebooks/teaching index 12c1a3711c..8ed7283af3 120000 --- a/notebooks/teaching +++ b/notebooks/teaching @@ -1 +1 @@ -../docs/source/teaching/ \ No newline at end of file +../docs/source/teaching/ diff --git a/notebooks/tutorials b/notebooks/tutorials index 8bfcbc5fe0..a499190667 120000 --- a/notebooks/tutorials +++ b/notebooks/tutorials @@ -1 +1 @@ -../docs/source/tutorials/ \ No newline at end of file +../docs/source/tutorials/ diff --git a/src/landlab/components/__init__.py b/src/landlab/components/__init__.py index 6a63465c45..cfb9358120 100644 --- a/src/landlab/components/__init__.py +++ b/src/landlab/components/__init__.py @@ -71,7 +71,7 @@ from .profiler import Profiler from .profiler import TrickleDownProfiler from .radiation import Radiation -from .river_flow_dynamics import river_flow_dynamics +from .river_flow_dynamics import RiverFlowDynamics from .sink_fill import SinkFiller from .sink_fill import SinkFillerBarnes from .soil_moisture import SoilInfiltrationGreenAmpt @@ -154,7 +154,7 @@ PrecipitationDistribution, Profiler, Radiation, - river_flow_dynamics, + RiverFlowDynamics, SedDepEroder, SedimentPulserAtLinks, SedimentPulserEachParcel, diff --git a/src/landlab/components/river_flow_dynamics/__init__.py b/src/landlab/components/river_flow_dynamics/__init__.py index ba4d1f1711..a93c0924d0 100644 --- a/src/landlab/components/river_flow_dynamics/__init__.py +++ b/src/landlab/components/river_flow_dynamics/__init__.py @@ -1,3 +1,3 @@ -from .river_flow_dynamics import river_flow_dynamics +from .river_flow_dynamics import RiverFlowDynamics -__all__ = ["river_flow_dynamics"] +__all__ = ["RiverFlowDynamics"] diff --git a/src/landlab/components/river_flow_dynamics/river_flow_dynamics.py b/src/landlab/components/river_flow_dynamics/river_flow_dynamics.py index d15b4a3d70..a8a7c6b4fd 100644 --- a/src/landlab/components/river_flow_dynamics/river_flow_dynamics.py +++ b/src/landlab/components/river_flow_dynamics/river_flow_dynamics.py @@ -6,104 +6,54 @@ Written by Sebastian Bernal and Angel Monsalve. -Last updated: October 18, 2023 - Examples -------- -River Flow Dynamics Simulation Example - -This example demonstrates how to simulate river flow dynamics using the Landlab library. - -First, import necessary libraries such as NumPy, Matplotlib, and Landlab components. +This example demonstrates basic usage of the RiverFlowDynamics component to simulate +a simple channel flow: >>> import numpy as np ->>> import matplotlib.pyplot as plt >>> from landlab import RasterModelGrid ->>> from landlab.components import river_flow_dynamics ->>> from landlab.io import read_esri_ascii ->>> from landlab.plot.imshow import imshow_grid - -Create a rectangular grid for flow dynamics calculations with specified dimensions -and cell size. - ->>> nRows = 20 ->>> nCols = 60 ->>> cellSize = 0.1 ->>> grid = RasterModelGrid((nRows, nCols), xy_spacing=(cellSize, cellSize)) - -Defining the Topography: Set up the initial topographic elevation for the grid, -creating a basic rectangular channel with a slope of 0.01. - ->>> te = grid.add_zeros("topographic__elevation", at="node") ->>> te += 0.059 - 0.01 * grid.x_of_node ->>> te[grid.y_of_node > 1.5] = 1.0 ->>> te[grid.y_of_node < 0.5] = 1.0 - -We could visualize the elevation profile using 'imshow_grid'. -imshow_grid(grid, "topographic__elevation") -Explore the grid's middle longitudinal section. - ->>> middleBedProfile = np.reshape(te, (nRows, nCols))[10, :] ->>> np.round(middleBedProfile, 3) -array([ 0.059, 0.058, 0.057, 0.056, 0.055, 0.054, 0.053, 0.052, - 0.051, 0.05 , 0.049, 0.048, 0.047, 0.046, 0.045, 0.044, - 0.043, 0.042, 0.041, 0.04 , 0.039, 0.038, 0.037, 0.036, - 0.035, 0.034, 0.033, 0.032, 0.031, 0.03 , 0.029, 0.028, - 0.027, 0.026, 0.025, 0.024, 0.023, 0.022, 0.021, 0.02 , - 0.019, 0.018, 0.017, 0.016, 0.015, 0.014, 0.013, 0.012, - 0.011, 0.01 , 0.009, 0.008, 0.007, 0.006, 0.005, 0.004, - 0.003, 0.002, 0.001, -0. ]) - -Instantiating the Component: Initialize the river_flow_dynamics component with -specified parameters, including the time step and Manning's roughness coefficient. -The grid will need some data to run the river_flow_dynamics component. -To check the names of the required inputs, use the 'input_var_names' class property. - ->>> river_flow_dynamics.input_var_names -('surface_water__depth', - 'surface_water__elevation', - 'surface_water__velocity', - 'topographic__elevation') +>>> from landlab.components import RiverFlowDynamics -To determine where these fields are mapped, use 'var_mapping'. +Create a small grid for demonstration purposes: ->>> river_flow_dynamics.var_mapping -(('surface_water__depth', 'node'), - ('surface_water__elevation', 'node'), - ('surface_water__velocity', 'link'), - ('topographic__elevation', 'node')) +>>> grid = RasterModelGrid((8, 6), xy_spacing=0.1) -Create fields of data for each of these input variables. The channel is initially empty, -so we create the water depth field. +Set up a sloped channel with elevated sides (slope of 0.01). ->>> h = grid.add_zeros("surface_water__depth", at="node") +>>> z = grid.add_zeros("topographic__elevation", at="node") +>>> z += 0.005 - 0.01 * grid.x_of_node +>>> z[grid.y_of_node > 0.5] = 1.0 +>>> z[grid.y_of_node < 0.2] = 1.0 -Water velocity is zero everywhere since there is no water yet. +Instantiating the Component. To check the names of the required inputs, use +the 'input_var_names' class property. ->>> vel = grid.add_zeros("surface_water__velocity", at="link") +>>> RiverFlowDynamics.input_var_names +('surface_water__depth', + 'surface_water__elevation', + 'surface_water__velocity', + 'topographic__elevation') -Calculate the initial water surface elevation from water depth and topographic elevation. +Initialize required fields: +>>> h = grid.add_zeros("surface_water__depth", at="node") +>>> vel = grid.add_zeros("surface_water__velocity", at="link") >>> wse = grid.add_zeros("surface_water__elevation", at="node") ->>> wse += h + te +>>> wse += h + z -Specify the nodes at which water enters the domain, and also the associated links. -These will serve as the inlet boundary conditions for water depth and velocity. -In this case, water flows from left to right at a depth of 0.5 meters with a velocity -of 0.45 m/s. +Set up inlet boundary conditions (left side of channel): +Water flows from left to right at a depth of 0.5 meters with a velocity of 0.45 m/s. ->>> fixed_entry_nodes = np.arange(300, 910, 60) +>>> fixed_entry_nodes = np.arange(12, 36, 6) >>> fixed_entry_links = grid.links_at_node[fixed_entry_nodes][:, 0] +>>> entry_nodes_h_values = np.full(4, 0.5) +>>> entry_links_vel_values = np.full(4, 0.45) -Set the fixed values for these entry nodes/links. - ->>> entry_nodes_h_values = np.full(11, 0.5) ->>> entry_links_vel_values = np.full(11, 0.45) - -Instantiate 'river_flow_dynamics' with the previously defined arguments. +Instantiate 'RiverFlowDynamics' ->>> rfd = river_flow_dynamics( +>>> rfd = RiverFlowDynamics( ... grid, ... dt=0.1, ... mannings_n=0.012, @@ -122,161 +72,17 @@ Examine the flow depth at the center of the channel after 10 seconds. -The expected flow depth is: - ->>> flow_depth_expected = np.array( -... [ -... 0.5, -... 0.491, -... 0.48, -... 0.473, -... 0.467, -... 0.464, -... 0.46, -... 0.458, -... 0.455, -... 0.454, -... 0.452, -... 0.45, -... 0.449, -... 0.448, -... 0.446, -... 0.445, -... 0.443, -... 0.442, -... 0.441, -... 0.439, -... 0.438, -... 0.437, -... 0.435, -... 0.434, -... 0.433, -... 0.431, -... 0.43, -... 0.428, -... 0.427, -... 0.425, -... 0.424, -... 0.422, -... 0.421, -... 0.419, -... 0.418, -... 0.416, -... 0.415, -... 0.413, -... 0.412, -... 0.41, -... 0.409, -... 0.407, -... 0.405, -... 0.404, -... 0.402, -... 0.401, -... 0.399, -... 0.397, -... 0.396, -... 0.394, -... 0.393, -... 0.391, -... 0.389, -... 0.388, -... 0.386, -... 0.384, -... 0.383, -... 0.381, -... 0.379, -... 0.378, -... ] -... ) - -The calculated flow_depth is: - ->>> flow_depth = np.reshape(grid["node"]["surface_water__depth"], (nRows, nCols))[10, :] - -The average (absolute) difference between predited and expected in percentage is: - ->>> np.round(np.abs(np.mean(flow_depth_expected - flow_depth)) * 100, 0) -0.0 +>>> flow_depth = np.reshape(grid["node"]["surface_water__depth"], (8, 6))[3, :] +>>> np.round(flow_depth, 3) +array([0.5 , 0.5 , 0.5 , 0.501, 0.502, 0.502]) And the velocity at links along the center of the channel. -The expected flow velocity is: - ->>> flow_velocity_expected = np.array( -... [ -... 0.45, -... 0.595, -... 0.694, -... 0.754, -... 0.795, -... 0.821, -... 0.838, -... 0.848, -... 0.855, -... 0.858, -... 0.86, -... 0.86, -... 0.858, -... 0.857, -... 0.858, -... 0.86, -... 0.864, -... 0.866, -... 0.866, -... 0.866, -... 0.866, -... 0.867, -... 0.869, -... 0.872, -... 0.874, -... 0.875, -... 0.876, -... 0.878, -... 0.88, -... 0.881, -... 0.882, -... 0.884, -... 0.885, -... 0.886, -... 0.888, -... 0.889, -... 0.89, -... 0.892, -... 0.893, -... 0.894, -... 0.896, -... 0.898, -... 0.898, -... 0.901, -... 0.901, -... 0.9, -... 0.904, -... 0.906, -... 0.902, -... 0.904, -... 0.91, -... 0.907, -... 0.904, -... 0.911, -... 0.911, -... 0.907, -... 0.909, -... 0.913, -... 0.914, -... ] -... ) - -The calculated flow velocity is: - ->>> linksAtCenter = grid.links_at_node[np.array(np.arange(600, 660))][:-1, 0] +>>> linksAtCenter = grid.links_at_node[np.array(np.arange(24, 30))][:-1, 0] >>> flow_velocity = grid["link"]["surface_water__velocity"][linksAtCenter] ->>> flow_velocity = np.round(flow_velocity, 3) - -The average (absolute) difference between predited and expected -flow velocity in percentage is: +>>> np.round(flow_velocity, 3) +array([0.45 , 0.457, 0.455, 0.452, 0.453]) ->>> np.round(np.abs(np.mean(flow_velocity_expected - flow_velocity)) * 100, 0) -0.0 """ import numpy as np @@ -285,7 +91,7 @@ from landlab import Component -class river_flow_dynamics(Component): +class RiverFlowDynamics(Component): """Simulate surface fluid flow based on Casulli and Cheng (1992). This Landlab component simulates surface fluid flow using the approximations of the @@ -306,7 +112,7 @@ class river_flow_dynamics(Component): https://doi.org/10.1002/fld.1650150602 """ - _name = "river_flow_dynamics" + _name = "RiverFlowDynamics" _unit_agnostic = False diff --git a/tests/components/river_flow_dynamics/test_river_flow_dynamics.py b/tests/components/river_flow_dynamics/test_river_flow_dynamics.py index 7a99079e7e..0c2459256f 100644 --- a/tests/components/river_flow_dynamics/test_river_flow_dynamics.py +++ b/tests/components/river_flow_dynamics/test_river_flow_dynamics.py @@ -1,14 +1,14 @@ """ -Unit tests for landlab.components.river_flow_dynamics.river_flow_dynamics +Unit tests for landlab.components.river_flow_dynamics.RiverFlowDynamics -last updated: 10/15/2023 +last updated: 11/11/2024 """ import numpy as np import pytest from landlab import RasterModelGrid -from landlab.components import river_flow_dynamics +from landlab.components import RiverFlowDynamics @pytest.fixture @@ -18,13 +18,13 @@ def rfd(): grid.add_zeros("surface_water__depth", at="node") grid.add_zeros("surface_water__elevation", at="node") grid.add_zeros("surface_water__velocity", at="link") - return river_flow_dynamics(grid) + return RiverFlowDynamics(grid) def test_name(rfd): """Test component name""" print(rfd.name) - assert rfd.name == "river_flow_dynamics" + assert rfd.name == "RiverFlowDynamics" def test_input_var_names(rfd): @@ -60,7 +60,7 @@ def test_initialization(): grid.add_zeros("surface_water__depth", at="node") grid.add_zeros("surface_water__elevation", at="node") grid.add_zeros("surface_water__velocity", at="link") - rfd = river_flow_dynamics(grid) + rfd = RiverFlowDynamics(grid) # Making sure fields have been created for field_name in rfd._info: @@ -71,7 +71,130 @@ def test_initialization(): # Re-initialize, this time with fields already existing in the grid # (this triggers the "if" instead of "else" in the field setup in init) - rfd = river_flow_dynamics(grid) + rfd = RiverFlowDynamics(grid) + + +def test_mass_conservation(): + """Test that water volume is conserved in a closed system.""" + grid = RasterModelGrid((10, 10), xy_spacing=0.1) + grid.add_zeros("topographic__elevation", at="node") + grid.add_zeros("surface_water__depth", at="node") + grid.add_zeros("surface_water__elevation", at="node") + grid.add_zeros("surface_water__velocity", at="link") + + # Create a bowl-shaped topography for more interesting flow + x = grid.x_of_node - 0.5 # Center at 0.5 + y = grid.y_of_node - 0.5 # Center at 0.5 + grid.at_node["topographic__elevation"] = 0.1 * (x**2 + y**2) + + # Set water surface elevation to 0.5 everywhere + grid.at_node["surface_water__elevation"][:] = 0.5 + + # Calculate water depth as difference between surface elevation and topography + topo = grid.at_node["topographic__elevation"] + water_surface = grid.at_node["surface_water__elevation"] + grid.at_node["surface_water__depth"] = np.maximum(water_surface - topo, 0.0) + + # Set up model with closed boundaries to ensure mass conservation + rfd = RiverFlowDynamics(grid) + + # Calculate initial volume + initial_volume = np.sum(grid.at_node["surface_water__depth"]) * grid.dx * grid.dy + + # Run model for several time steps + n_steps = 100 + volumes = np.zeros(n_steps + 1) + volumes[0] = initial_volume + + for i in range(n_steps): + rfd.run_one_step() + volumes[i + 1] = ( + np.sum(grid.at_node["surface_water__depth"]) * grid.dx * grid.dy + ) + + # Calculate volume change + volume_changes = np.abs(volumes - initial_volume) / initial_volume + max_volume_change = np.max(volume_changes) + + # Assert volume is conserved within a more realistic tolerance (0.1%) + assert max_volume_change < 1e-3 + + # Check for physically reasonable results + depths = grid.at_node["surface_water__depth"] + initial_max_depth = np.max(grid.at_node["surface_water__depth"]) + assert np.all(depths >= 0) + assert np.max(depths) < initial_max_depth * 1.001 + + +def test_open_boundaries(): + """Test that water properly exits through open boundaries.""" + grid = RasterModelGrid((10, 10), xy_spacing=0.1) + grid.add_zeros("topographic__elevation", at="node") + grid.add_zeros("surface_water__depth", at="node") + grid.add_zeros("surface_water__elevation", at="node") + grid.add_zeros("surface_water__velocity", at="link") + + # Sloped channel + grid.at_node["topographic__elevation"] = grid.x_of_node * 0.05 + + # Set up inflow conditions + fixed_entry_nodes = grid.nodes_at_left_edge + fixed_entry_links = grid.links_at_node[fixed_entry_nodes][:, 0] + entry_nodes_h_values = 0.5 * np.ones_like(fixed_entry_nodes) + entry_links_vel_values = 0.3 * np.ones_like(fixed_entry_links) + + rfd = RiverFlowDynamics( + grid, + fixed_entry_nodes=fixed_entry_nodes, + fixed_entry_links=fixed_entry_links, + entry_nodes_h_values=entry_nodes_h_values, + entry_links_vel_values=entry_links_vel_values, + ) + + # Run to steady state + for _ in range(200): + rfd.run_one_step() + + # Check that water exits smoothly (no backing up at boundary) + right_depths = grid.at_node["surface_water__depth"][grid.nodes_at_right_edge] + near_right_depths = grid.at_node["surface_water__depth"][ + grid.nodes_at_right_edge - 1 + ] + + # Water depth should decrease or stay similar near boundary + assert np.all(right_depths <= near_right_depths * 1.1) + + +def test_still_water(): + """Test that still water stays still.""" + grid = RasterModelGrid((10, 10), xy_spacing=0.1) + grid.add_zeros("topographic__elevation", at="node") + grid.add_zeros("surface_water__depth", at="node") + grid.add_zeros("surface_water__elevation", at="node") + grid.add_zeros("surface_water__velocity", at="link") + + # Flat bottom with uniform water depth + grid.at_node["surface_water__depth"][:] = 0.5 + grid.at_node["surface_water__elevation"][:] = 0.5 + + rfd = RiverFlowDynamics(grid) + + # Initial conditions + initial_volume = np.sum(grid.at_node["surface_water__depth"]) * grid.dx * grid.dy + + # Run model + rfd.run_one_step() + + # Check volume conservation + final_volume = np.sum(grid.at_node["surface_water__depth"]) * grid.dx * grid.dy + volume_change = abs(final_volume - initial_volume) / initial_volume + assert volume_change < 1e-4 # Allow 0.01% volume change + + # Check surface stays flat (no gradients) + depth_gradients = np.abs( + np.diff(grid.at_node["surface_water__depth"].reshape(10, 10), axis=1) + ) + assert np.all(depth_gradients < 1e-4) # Maximum gradient of 0.1 mm per cell def test_run_one_step(): @@ -92,7 +215,7 @@ def test_run_one_step(): # We specify the time step duration and we run it dt = 0.1 - rfd = river_flow_dynamics( + rfd = RiverFlowDynamics( grid, dt=dt, fixed_entry_nodes=fixed_entry_nodes, @@ -129,3 +252,257 @@ def test_run_one_step(): np.testing.assert_array_almost_equal( water_depth_solution, water_depth_obtained, decimal=3 ) + + +def test_downhill_flow(): + """Test that water flows downhill.""" + grid = RasterModelGrid((10, 10), xy_spacing=0.1) + grid.add_zeros("topographic__elevation", at="node") + grid.add_zeros("surface_water__depth", at="node") + grid.add_zeros("surface_water__elevation", at="node") + grid.add_zeros("surface_water__velocity", at="link") + + # Create sloped surface + grid.at_node["topographic__elevation"] = grid.x_of_node * 0.1 + + # Set water surface elevation at upstream end + upstream_water_level = 0.5 + grid.at_node["surface_water__elevation"][ + grid.nodes_at_left_edge + ] = upstream_water_level + + # Calculate initial water depth at upstream end + upstream_topo = grid.at_node["topographic__elevation"][grid.nodes_at_left_edge] + grid.at_node["surface_water__depth"][grid.nodes_at_left_edge] = ( + upstream_water_level - upstream_topo + ) + + # Set up inflow conditions + fixed_entry_nodes = grid.nodes_at_left_edge + fixed_entry_links = grid.links_at_node[fixed_entry_nodes][:, 0] + entry_links_vel_values = 0.3 * np.ones_like(fixed_entry_links) + + rfd = RiverFlowDynamics( + grid, + fixed_entry_nodes=fixed_entry_nodes, + fixed_entry_links=fixed_entry_links, + entry_nodes_h_values=grid.at_node["surface_water__depth"][fixed_entry_nodes], + entry_links_vel_values=entry_links_vel_values, + ) + + # Run model + for _ in range(50): + rfd.run_one_step() + + # Check mean flow direction + horizontal_links = grid.horizontal_links + mean_velocity = np.mean(grid.at_link["surface_water__velocity"][horizontal_links]) + assert mean_velocity > 0 # Overall flow should be positive (downhill) + + # Check water depths are physically reasonable + assert np.all(grid.at_node["surface_water__depth"] >= 0) + + # Check water surface elevation is consistent with depth and topography + # but allow for small numerical differences + calculated_elevation = ( + grid.at_node["surface_water__depth"] + grid.at_node["topographic__elevation"] + ) + elevation_difference = np.abs( + grid.at_node["surface_water__elevation"] - calculated_elevation + ) + assert np.max(elevation_difference) < 0.1 # Allow differences up to 10cm + + # Check that water has moved downstream + outlet_depths = grid.at_node["surface_water__depth"][grid.nodes_at_right_edge] + assert np.any(outlet_depths > 0.01) + + +def test_time_step_sensitivity(): + """Test that results are not heavily dependent on time step choice for + relatively similar systems.""" + grid1 = RasterModelGrid((10, 10), xy_spacing=0.1) + grid2 = RasterModelGrid((10, 10), xy_spacing=0.1) + + # Set up identical grids + for grid in [grid1, grid2]: + grid.add_zeros("topographic__elevation", at="node") + grid.add_zeros("surface_water__depth", at="node") + grid.add_zeros("surface_water__elevation", at="node") + grid.add_zeros("surface_water__velocity", at="link") + + # Create a sloped surface + grid.at_node["topographic__elevation"] = grid.x_of_node * 0.1 + + # Set the water surface elevation + water_level = 0.5 + grid.at_node["surface_water__elevation"][grid.nodes_at_left_edge] = water_level + + # Calculate initial water depth + topo = grid.at_node["topographic__elevation"][grid.nodes_at_left_edge] + grid.at_node["surface_water__depth"][grid.nodes_at_left_edge] = ( + water_level - topo + ) + + # Set up identical boundary conditions for both grids + fixed_entry_nodes = grid1.nodes_at_left_edge + fixed_entry_links = grid1.links_at_node[fixed_entry_nodes][:, 0] + entry_nodes_h_values = grid1.at_node["surface_water__depth"][fixed_entry_nodes] + entry_links_vel_values = 0.3 * np.ones_like(fixed_entry_links) + + # Create two models with different time steps + rfd1 = RiverFlowDynamics( + grid1, + dt=0.1, + fixed_entry_nodes=fixed_entry_nodes, + fixed_entry_links=fixed_entry_links, + entry_nodes_h_values=entry_nodes_h_values, + entry_links_vel_values=entry_links_vel_values, + ) + + rfd2 = RiverFlowDynamics( + grid2, + dt=0.05, + fixed_entry_nodes=fixed_entry_nodes, + fixed_entry_links=fixed_entry_links, + entry_nodes_h_values=entry_nodes_h_values, + entry_links_vel_values=entry_links_vel_values, + ) + + # Run to same total time + for _ in range(100): + rfd1.run_one_step() + for _ in range(200): + rfd2.run_one_step() + + # Compare results with reasonable tolerances + depth_difference = np.abs( + grid1.at_node["surface_water__depth"] - grid2.at_node["surface_water__depth"] + ) + assert np.max(depth_difference) < 0.1 # Max difference in depths + + # Check both solutions maintain physical consistency + for grid in [grid1, grid2]: + # Water depth should be non-negative + assert np.all(grid.at_node["surface_water__depth"] >= 0) + + # Check water surface elevation consistency + calculated_elevation = ( + grid.at_node["surface_water__depth"] + + grid.at_node["topographic__elevation"] + ) + elevation_difference = np.abs( + grid.at_node["surface_water__elevation"] - calculated_elevation + ) + assert np.max(elevation_difference) < 0.1 + + # Check water has propagated downstream + outlet_depths = grid.at_node["surface_water__depth"][grid.nodes_at_right_edge] + assert np.any(outlet_depths > 0.01) + + +def test_numerical_stability(): + """Test numerical stability under various conditions.""" + + def create_test_grid(): + """Helper function to create a grid with initial conditions.""" + grid = RasterModelGrid((10, 10), xy_spacing=0.1) + grid.add_zeros("topographic__elevation", at="node") + grid.add_zeros("surface_water__depth", at="node") + grid.add_zeros("surface_water__elevation", at="node") + grid.add_zeros("surface_water__velocity", at="link") + + # Create challenging conditions + grid.at_node["topographic__elevation"] = grid.x_of_node * 0.1 + + # Set initial conditions + grid.at_node["surface_water__depth"][:] = 0.5 + grid.at_node["surface_water__elevation"] = ( + grid.at_node["surface_water__depth"] + + grid.at_node["topographic__elevation"] + ) + + return grid + + time_steps = [0.1, 0.05, 0.01] + max_velocities = [] + final_depths = [] + + for dt in time_steps: + grid_test = create_test_grid() + rfd = RiverFlowDynamics(grid_test, dt=dt) + + try: + for _ in range(int(5.0 / dt)): + rfd.run_one_step() + + assert np.all(np.isfinite(grid_test.at_node["surface_water__depth"])) + assert np.all(grid_test.at_node["surface_water__depth"] >= 0) + assert np.all(np.isfinite(grid_test.at_link["surface_water__velocity"])) + + max_velocities.append( + np.max(np.abs(grid_test.at_link["surface_water__velocity"])) + ) + final_depths.append(grid_test.at_node["surface_water__depth"].copy()) + + except AssertionError: + raise AssertionError(f"Instability detected with dt={dt}") + + depth_variations = [ + np.max(np.abs(d1 - d2)) for d1, d2 in zip(final_depths[:-1], final_depths[1:]) + ] + assert np.all(np.array(depth_variations) < 0.1) + assert np.all(np.array(max_velocities) < 1.0) + + +def test_analytical_solution(): + """Test numerical and analytical solutions.""" + + # Creating a grid + grid = RasterModelGrid((6, 6), xy_spacing=0.1) + grid.add_zeros("topographic__elevation", at="node") + grid.add_zeros("surface_water__depth", at="node") + grid.add_zeros("surface_water__elevation", at="node") + grid.add_zeros("surface_water__velocity", at="link") + + # We set fixed boundary conditions, specifying the nodes and links in which + # the water is flowing into the grid + fixed_entry_nodes = grid.nodes_at_left_edge + fixed_entry_links = grid.links_at_node[fixed_entry_nodes][:, 0] + + # We set the fixed values in the entry nodes/links + entry_nodes_h_values = 0.5 * np.ones_like(fixed_entry_nodes) + # entry_nodes_h_values = np.repeat(0.5, 6) + entry_links_vel_values = 0.6 * np.ones_like(fixed_entry_links) + # entry_links_vel_values = np.repeat(0.6, 6) + + # We specify the time step duration and we run it + dt = 0.1 + rfd = RiverFlowDynamics( + grid, + dt=dt, + mannings_n=0.0010, + fixed_entry_nodes=fixed_entry_nodes, + fixed_entry_links=fixed_entry_links, + entry_nodes_h_values=entry_nodes_h_values, + entry_links_vel_values=entry_links_vel_values, + ) + + for _ in range(300): + rfd.run_one_step() + + # Comparing the uniform outlet with the obtained water depth + water_depth_solution = 0.5 * np.ones_like(fixed_entry_nodes) + water_depth_obtained = grid.at_node["surface_water__depth"][ + grid.nodes_at_right_edge + ] + + water_depth_difference = np.abs(water_depth_solution - water_depth_obtained) + assert np.max(water_depth_difference) < 0.001 # Max difference in vel + + vel_solution = 0.6 * np.ones_like(fixed_entry_links)[1:5] + vel_obtained = grid.at_link["surface_water__velocity"][ + grid.links_at_node[grid.nodes_at_right_edge][:, 2] + ][1:5] + vel_difference = np.abs(vel_solution - vel_obtained) + + assert np.max(vel_difference) < 0.01 # Max difference in vel