diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 3172ea84..b8d0db84 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,6 +1,6 @@ repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.5.0 + rev: v5.0.0 hooks: - id: trailing-whitespace - id: end-of-file-fixer @@ -8,17 +8,21 @@ repos: - id: debug-statements - id: mixed-line-ending - repo: https://github.com/PyCQA/isort - rev: 5.12.0 + rev: 5.13.2 hooks: - id: isort args: ["--profile", "black"] - repo: https://github.com/psf/black - rev: 23.9.1 + rev: 24.10.0 hooks: - id: black args: [--safe, --quiet] - repo: https://github.com/pre-commit/mirrors-clang-format - rev: v17.0.2 + rev: v19.1.1 hooks: - id: clang-format args: [--style=file] + +ci: + autofix_prs: false + autoupdate_schedule: monthly diff --git a/doc/source/examples/catchment_py.ipynb b/doc/source/examples/catchment_py.ipynb index a4471387..a5ba1e4f 100644 --- a/doc/source/examples/catchment_py.ipynb +++ b/doc/source/examples/catchment_py.ipynb @@ -1,319 +1,270 @@ { "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Catchment (TriMesh) [Py]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "A simple example simulating the evolution of a single synthetic catchment on a triangular (irregular) mesh.\n", - "\n", - "The local rate of elevation change, {math}`\\partial h/\\partial t`, is governed by bedrock river channel erosion (stream-power law).\n", - "\n", - "```{math}\n", - ":label: eq_catchment\n", - "\\frac{\\partial h}{\\partial t} = - K_f A^m (\\nabla h)^n\n", - "```\n", - "\n", - "where {math}`A` is the drainage area (i.e., the total upslope area from which the water flow converges) and {math}`K_f`, {math}`m` and {math}`n` are parameters (uniform in space and time).\n", - "\n", - "The initial topography is a nearly flat plateau of a given elevation. The base level (catchment outlet) is lowered during the simulation at a constant rate." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import matplotlib.tri as mpltri\n", - "import pygalmesh\n", - "\n", - "import fastscapelib as fs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "remove-cell" - ] - }, - "outputs": [], - "source": [ - "# Theme that looks reasonably fine on both dark/light modes\n", - "import matplotlib\n", - "matplotlib.style.use('Solarize_Light2')\n", - "matplotlib.rcParams['axes.grid'] = False" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup the Mesh, Flow Graph and Eroders\n", - "\n", - "Import the coordinates of the catchment boundaries from a file and create a new irregular mesh (with quasi-uniform resolution) using a constrained triangulation algorithm available in the [pygalmesh](https://github.com/meshpro/pygalmesh) library." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "basin_poly = np.load(\"basin.npy\")\n", - "n_poly = basin_poly.shape[0]\n", - "constraints = list([list(s) for s in zip(range(n_poly), range(1, n_poly))])\n", - "constraints[-1][1] = 0\n", - "\n", - "approx_resolution = 300.0\n", - "\n", - "mesh = pygalmesh.generate_2d(\n", - " basin_poly,\n", - " constraints,\n", - " max_edge_size=approx_resolution,\n", - " num_lloyd_steps=2,\n", - ")\n", - "\n", - "points = mesh.points\n", - "triangles = mesh.cells[0].data" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Determine the base level node, which corresponds to the catchment outlet (a boundary node)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "outlet_idx = np.argmax(points[:, 1])\n", - "base_levels = {outlet_idx: fs.NodeStatus.FIXED_VALUE}" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create a new {py:class}`~fastscapelib.TriMesh` object and pass the existing mesh data." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "grid = fs.TriMesh(points, triangles, base_levels)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create a {py:class}`~fastscapelib.FlowGraph` object with single direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n", - "\n", - "By default, base level nodes are set from fixed value boundary conditions (the outlet node set above in this example)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setup the bedrock channel eroder with a given set of parameter values." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "spl_eroder = fs.SPLEroder(\n", - " flow_graph,\n", - " k_coef=2e-4,\n", - " area_exp=0.4,\n", - " slope_exp=1,\n", - " tolerance=1e-5,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup Initial Conditions and External Forcing\n", - "\n", - "Create a flat (+ random perturbations) elevated surface topography as initial conditions. Also initialize the array for drainage area." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "init_elevation = np.random.uniform(0, 1, size=grid.shape) + 500.0\n", - "elevation = init_elevation\n", - "drainage_area = np.empty_like(elevation)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The base level lowering rate:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "lowering_rate = 1e-4" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run the Model\n", - "\n", - "Run the model for one hundred time steps." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dt = 2e4\n", - "nsteps = 100\n", - "\n", - "for step in range(100):\n", - " # base level (outlet) lowering \n", - " elevation[outlet_idx] -= lowering_rate * dt\n", - " \n", - " # flow routing\n", - " flow_graph.update_routes(elevation)\n", - " \n", - " # flow accumulation (drainage area)\n", - " flow_graph.accumulate(drainage_area, 1.0)\n", - " \n", - " # channel erosion (SPL)\n", - " spl_erosion = spl_eroder.erode(elevation, drainage_area, dt)\n", - " \n", - " # update topography\n", - " elevation -= spl_erosion" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Plot Outputs and Other Diagnostics" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x, y = points.transpose()\n", - "plot_tri = mpltri.Triangulation(x, y, triangles)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Topographic elevation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(13, 13))\n", - "ax.set_aspect('equal')\n", - "cf = ax.tripcolor(plot_tri, elevation)\n", - "plt.colorbar(cf)\n", - "ax.triplot(plot_tri, linewidth=0.2, c=\"k\");" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Drainage area (log)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(13, 13))\n", - "ax.set_aspect('equal')\n", - "cf = ax.tricontourf(plot_tri, np.log(drainage_area), cmap=plt.cm.Blues)\n", - "plt.colorbar(cf);\n", - "ax.triplot(plot_tri, linewidth=0.2, c=\"k\");" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python [conda env:fastscapelib_dev]", - "language": "python", - "name": "conda-env-fastscapelib_dev-py" - }, - "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.9.10" - }, - "vscode": { - "interpreter": { - "hash": "828c8bd590fe165438ea9436ea3075ad1e985d4b674461d5f3ca65e5df0d3dd5" - } - } - }, - "nbformat": 4, - "nbformat_minor": 4 + { "cell_type": "markdown", "metadata": {}, "source": ["# Catchment (TriMesh) [Py]"] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A simple example simulating the evolution of a single synthetic catchment on a triangular (irregular) mesh.\n", + "\n", + "The local rate of elevation change, {math}`\\partial h/\\partial t`, is governed by bedrock river channel erosion (stream-power law).\n", + "\n", + "```{math}\n", + ":label: eq_catchment\n", + "\\frac{\\partial h}{\\partial t} = - K_f A^m (\\nabla h)^n\n", + "```\n", + "\n", + "where {math}`A` is the drainage area (i.e., the total upslope area from which the water flow converges) and {math}`K_f`, {math}`m` and {math}`n` are parameters (uniform in space and time).\n", + "\n", + "The initial topography is a nearly flat plateau of a given elevation. The base level (catchment outlet) is lowered during the simulation at a constant rate." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.tri as mpltri\n", + "import pygalmesh\n", + "\n", + "import fastscapelib as fs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { "tags": ["remove-cell"] }, + "outputs": [], + "source": [ + "# Theme that looks reasonably fine on both dark/light modes\n", + "import matplotlib\n", + "matplotlib.style.use('Solarize_Light2')\n", + "matplotlib.rcParams['axes.grid'] = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup the Mesh, Flow Graph and Eroders\n", + "\n", + "Import the coordinates of the catchment boundaries from a file and create a new irregular mesh (with quasi-uniform resolution) using a constrained triangulation algorithm available in the [pygalmesh](https://github.com/meshpro/pygalmesh) library." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "basin_poly = np.load(\"basin.npy\")\n", + "n_poly = basin_poly.shape[0]\n", + "constraints = list([list(s) for s in zip(range(n_poly), range(1, n_poly))])\n", + "constraints[-1][1] = 0\n", + "\n", + "approx_resolution = 300.0\n", + "\n", + "mesh = pygalmesh.generate_2d(\n", + " basin_poly,\n", + " constraints,\n", + " max_edge_size=approx_resolution,\n", + " num_lloyd_steps=2,\n", + ")\n", + "\n", + "points = mesh.points\n", + "triangles = mesh.cells[0].data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Determine the base level node, which corresponds to the catchment outlet (a boundary node)."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "outlet_idx = np.argmax(points[:, 1])\n", + "base_levels = {outlet_idx: fs.NodeStatus.FIXED_VALUE}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Create a new {py:class}`~fastscapelib.TriMesh` object and pass the existing mesh data."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": ["grid = fs.TriMesh(points, triangles, base_levels)"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a {py:class}`~fastscapelib.FlowGraph` object with single direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n", + "\n", + "By default, base level nodes are set from fixed value boundary conditions (the outlet node set above in this example)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": + ["flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["Setup the bedrock channel eroder with a given set of parameter values."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spl_eroder = fs.SPLEroder(\n", + " flow_graph,\n", + " k_coef=2e-4,\n", + " area_exp=0.4,\n", + " slope_exp=1,\n", + " tolerance=1e-5,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup Initial Conditions and External Forcing\n", + "\n", + "Create a flat (+ random perturbations) elevated surface topography as initial conditions. Also initialize the array for drainage area." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "init_elevation = np.random.uniform(0, 1, size=grid.shape) + 500.0\n", + "elevation = init_elevation\n", + "drainage_area = np.empty_like(elevation)" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["The base level lowering rate:"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": ["lowering_rate = 1e-4"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Run the Model\n", "\n", "Run the model for one hundred time steps."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = 2e4\n", + "nsteps = 100\n", + "\n", + "for step in range(100):\n", + " # base level (outlet) lowering \n", + " elevation[outlet_idx] -= lowering_rate * dt\n", + " \n", + " # flow routing\n", + " flow_graph.update_routes(elevation)\n", + " \n", + " # flow accumulation (drainage area)\n", + " flow_graph.accumulate(drainage_area, 1.0)\n", + " \n", + " # channel erosion (SPL)\n", + " spl_erosion = spl_eroder.erode(elevation, drainage_area, dt)\n", + " \n", + " # update topography\n", + " elevation -= spl_erosion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Plot Outputs and Other Diagnostics"] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": + ["x, y = points.transpose()\n", "plot_tri = mpltri.Triangulation(x, y, triangles)"] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Topographic elevation"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(13, 13))\n", + "ax.set_aspect('equal')\n", + "cf = ax.tripcolor(plot_tri, elevation)\n", + "plt.colorbar(cf)\n", + "ax.triplot(plot_tri, linewidth=0.2, c=\"k\");" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Drainage area (log)"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(13, 13))\n", + "ax.set_aspect('equal')\n", + "cf = ax.tricontourf(plot_tri, np.log(drainage_area), cmap=plt.cm.Blues)\n", + "plt.colorbar(cf);\n", + "ax.triplot(plot_tri, linewidth=0.2, c=\"k\");" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:fastscapelib_dev]", + "language": "python", + "name": "conda-env-fastscapelib_dev-py" + }, + "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.9.10" + }, + "vscode": { + "interpreter": + { "hash": "828c8bd590fe165438ea9436ea3075ad1e985d4b674461d5f3ca65e5df0d3dd5" } + } + }, + "nbformat": 4, + "nbformat_minor": 4 } diff --git a/doc/source/examples/escarpment_py.ipynb b/doc/source/examples/escarpment_py.ipynb index 033d9520..5004d9c5 100644 --- a/doc/source/examples/escarpment_py.ipynb +++ b/doc/source/examples/escarpment_py.ipynb @@ -1,252 +1,217 @@ { "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Escarpment [Py]\n", - "\n", - "A simple example simulating the evolution of an escarpment (on a 2D raster grid) under the action of bedrock channel and hillslope erosion.\n", - "\n", - "This example is pretty similar to {doc}`mountain_py`, although solving equation {eq}`eq_mountain` without the uplift term {math}`U` on a semi-infinite (periodic) domain. It also starts from a different initial topograhic surface and routes flow using a multiple direction approach." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import fastscapelib as fs\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "remove-cell" - ] - }, - "outputs": [], - "source": [ - "# Theme that looks reasonably fine on both dark/light modes\n", - "import matplotlib\n", - "matplotlib.style.use('Solarize_Light2')\n", - "matplotlib.rcParams['axes.grid'] = False" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup the Grid, Flow Graph and Eroders\n", - "\n", - "Create a {py:class}`~fastscapelib.RasterGrid` of 101x201 nodes with a total length of 10 km in y (rows) and 20 km in x (columns).\n", - "\n", - "Set fixed-value boundary at the left border, free (core) boundary at the right border and reflective boundaries for the top-down borders." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bs = [\n", - " fs.NodeStatus.FIXED_VALUE,\n", - " fs.NodeStatus.CORE,\n", - " fs.NodeStatus.LOOPED,\n", - " fs.NodeStatus.LOOPED,\n", - "]\n", - "\n", - "grid = fs.RasterGrid.from_length([101, 201], [1e4, 2e4], bs)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create a {py:class}`~fastscapelib.FlowGraph` object with multiple direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n", - "\n", - "By default, base level nodes are set from fixed value boundary conditions (all nodes on the left border in this example)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "flow_graph = fs.FlowGraph(\n", - " grid,\n", - " [fs.SingleFlowRouter(), fs.MSTSinkResolver(), fs.MultiFlowRouter(1.1)],\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setup eroder classes (bedrock channel + hillslope) with a given set of parameter values." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "spl_eroder = fs.SPLEroder(\n", - " flow_graph,\n", - " k_coef=1e-4,\n", - " area_exp=0.4,\n", - " slope_exp=1,\n", - " tolerance=1e-5,\n", - ")\n", - "\n", - "diffusion_eroder = fs.DiffusionADIEroder(grid, 0.01)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup Initial Conditions and External Forcing" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setup initial topography as two nearly flat plateaus (+ small random perturbations) separated by a steep escarpment. Also initialize the array for drainage area." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "rng = np.random.Generator(np.random.PCG64(1234))\n", - "\n", - "init_elevation = rng.uniform(0, 1, size=grid.shape)\n", - "init_elevation[:, 100:] += 400\n", - "\n", - "elevation = init_elevation\n", - "drainage_area = np.empty_like(elevation)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run the Model" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Run the model for a few dozens of time steps (total simulation time: 100k years)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dt = 2e3\n", - "nsteps = 50\n", - "\n", - "for step in range(nsteps):\n", - " # flow routing\n", - " flow_graph.update_routes(elevation)\n", - " \n", - " # flow accumulation (drainage area)\n", - " flow_graph.accumulate(drainage_area, 1.0)\n", - " \n", - " # apply channel erosion then hillslope diffusion\n", - " spl_erosion = spl_eroder.erode(elevation, drainage_area, dt)\n", - " diff_erosion = diffusion_eroder.erode(elevation - spl_erosion, dt)\n", - " \n", - " # update topography\n", - " elevation = elevation - spl_erosion - diff_erosion\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Plot Outputs and Other Diagnostics\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Topographic elevation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(12, 6))\n", - "plt.imshow(elevation)\n", - "plt.colorbar();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Drainage area (log)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(12, 6))\n", - "plt.imshow(np.log(drainage_area), cmap=plt.cm.Blues)\n", - "plt.colorbar();" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "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.9.10" - } - }, - "nbformat": 4, - "nbformat_minor": 4 + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Escarpment [Py]\n", + "\n", + "A simple example simulating the evolution of an escarpment (on a 2D raster grid) under the action of bedrock channel and hillslope erosion.\n", + "\n", + "This example is pretty similar to {doc}`mountain_py`, although solving equation {eq}`eq_mountain` without the uplift term {math}`U` on a semi-infinite (periodic) domain. It also starts from a different initial topograhic surface and routes flow using a multiple direction approach." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import fastscapelib as fs\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { "tags": ["remove-cell"] }, + "outputs": [], + "source": [ + "# Theme that looks reasonably fine on both dark/light modes\n", + "import matplotlib\n", + "matplotlib.style.use('Solarize_Light2')\n", + "matplotlib.rcParams['axes.grid'] = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup the Grid, Flow Graph and Eroders\n", + "\n", + "Create a {py:class}`~fastscapelib.RasterGrid` of 101x201 nodes with a total length of 10 km in y (rows) and 20 km in x (columns).\n", + "\n", + "Set fixed-value boundary at the left border, free (core) boundary at the right border and reflective boundaries for the top-down borders." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bs = [\n", + " fs.NodeStatus.FIXED_VALUE,\n", + " fs.NodeStatus.CORE,\n", + " fs.NodeStatus.LOOPED,\n", + " fs.NodeStatus.LOOPED,\n", + "]\n", + "\n", + "grid = fs.RasterGrid.from_length([101, 201], [1e4, 2e4], bs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a {py:class}`~fastscapelib.FlowGraph` object with multiple direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n", + "\n", + "By default, base level nodes are set from fixed value boundary conditions (all nodes on the left border in this example)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "flow_graph = fs.FlowGraph(\n", + " grid,\n", + " [fs.SingleFlowRouter(), fs.MSTSinkResolver(), fs.MultiFlowRouter(1.1)],\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Setup eroder classes (bedrock channel + hillslope) with a given set of parameter values."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spl_eroder = fs.SPLEroder(\n", + " flow_graph,\n", + " k_coef=1e-4,\n", + " area_exp=0.4,\n", + " slope_exp=1,\n", + " tolerance=1e-5,\n", + ")\n", + "\n", + "diffusion_eroder = fs.DiffusionADIEroder(grid, 0.01)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Setup Initial Conditions and External Forcing"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Setup initial topography as two nearly flat plateaus (+ small random perturbations) separated by a steep escarpment. Also initialize the array for drainage area."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rng = np.random.Generator(np.random.PCG64(1234))\n", + "\n", + "init_elevation = rng.uniform(0, 1, size=grid.shape)\n", + "init_elevation[:, 100:] += 400\n", + "\n", + "elevation = init_elevation\n", + "drainage_area = np.empty_like(elevation)" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["## Run the Model"] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Run the model for a few dozens of time steps (total simulation time: 100k years)."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = 2e3\n", + "nsteps = 50\n", + "\n", + "for step in range(nsteps):\n", + " # flow routing\n", + " flow_graph.update_routes(elevation)\n", + " \n", + " # flow accumulation (drainage area)\n", + " flow_graph.accumulate(drainage_area, 1.0)\n", + " \n", + " # apply channel erosion then hillslope diffusion\n", + " spl_erosion = spl_eroder.erode(elevation, drainage_area, dt)\n", + " diff_erosion = diffusion_eroder.erode(elevation - spl_erosion, dt)\n", + " \n", + " # update topography\n", + " elevation = elevation - spl_erosion - diff_erosion\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Plot Outputs and Other Diagnostics\n"] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Topographic elevation"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 6))\n", + "plt.imshow(elevation)\n", + "plt.colorbar();" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Drainage area (log)"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 6))\n", + "plt.imshow(np.log(drainage_area), cmap=plt.cm.Blues)\n", + "plt.colorbar();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.9.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 } diff --git a/doc/source/examples/inner_base_levels_py.ipynb b/doc/source/examples/inner_base_levels_py.ipynb index b4314a15..5dc77bb0 100644 --- a/doc/source/examples/inner_base_levels_py.ipynb +++ b/doc/source/examples/inner_base_levels_py.ipynb @@ -1,299 +1,256 @@ { "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Inner Base Levels [Py]\n", - "\n", - "An example very similar to {doc}`mountain_py` but setting the base level nodes inside the raster grid instead of on the boundaries. The modeled domain is thus infinite (fully periodic), which is usefull for global (planetary) scale simulations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from random import random, uniform\n", - "\n", - "import fastscapelib as fs\n", - "import numpy as np\n", - "import matplotlib\n", - "import matplotlib.pyplot as plt" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "remove-cell" - ] - }, - "outputs": [], - "source": [ - "# Theme that looks reasonably fine on both dark/light modes\n", - "matplotlib.style.use('Solarize_Light2')\n", - "matplotlib.rcParams['axes.grid'] = False" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup the Grid, Flow Graph and Eroders\n", - "\n", - "Create a {py:class}`~fastscapelib.RasterGrid` of 201x201 nodes with a total length of 50 km in both y (rows) and x (columns).\n", - "\n", - "Set looped (reflective) boundary conditions at all border nodes. Also set fixed-value \"boundary\" for a given number of base level nodes randomly selected inside the domain." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "n_base_levels = 20\n", - "base_level_row = np.random.uniform(1, 200, n_base_levels).astype(\"int\")\n", - "base_level_col = np.random.uniform(1, 200, n_base_levels).astype(\"int\")\n", - "\n", - "base_levels = {\n", - " (i, j): fs.NodeStatus.FIXED_VALUE\n", - " for i, j in zip(base_level_row, base_level_col)\n", - "}\n", - "\n", - "bs = fs.NodeStatus.LOOPED\n", - "grid = fs.RasterGrid.from_length([201, 201], [5e4, 5e4], bs, base_levels)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create a {py:class}`~fastscapelib.FlowGraph` object with single direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n", - "\n", - "By default, base level nodes are set from fixed value boundary conditions (random inner nodes in this example)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setup eroder classes (bedrock channel + hillslope) with a given set of parameter values." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "spl_eroder = fs.SPLEroder(\n", - " flow_graph,\n", - " k_coef=2e-4,\n", - " area_exp=0.4,\n", - " slope_exp=1,\n", - " tolerance=1e-5,\n", - ")\n", - "\n", - "diffusion_eroder = fs.DiffusionADIEroder(grid, 0.01)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup Initial Conditions and External Forcing" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create a flat (+ random perturbations) surface topography as initial conditions. Also initialize the array for drainage area." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "rng = np.random.Generator(np.random.PCG64(1234))\n", - "\n", - "init_elevation = rng.uniform(0, 1, size=grid.shape)\n", - "\n", - "elevation = init_elevation\n", - "drainage_area = np.empty_like(elevation)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Set upflit rate as uniform (fixed value) within the domain and to zero at all grid boundaries." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "uplift_rate = np.where(\n", - " grid.nodes_status() == fs.NodeStatus.FIXED_VALUE.value, 0, 1e-3\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run the Model" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Run the model for a few dozens of time steps (total simulation time: 1M years)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dt = 2e4\n", - "nsteps = 50\n", - "\n", - "for step in range(nsteps):\n", - " # uplift (no uplift at fixed elevation boundaries)\n", - " uplifted_elevation = elevation + dt * uplift_rate\n", - " \n", - " # flow routing\n", - " filled_elevation = flow_graph.update_routes(uplifted_elevation)\n", - " \n", - " # flow accumulation (drainage area)\n", - " flow_graph.accumulate(drainage_area, 1.0)\n", - " \n", - " # apply channel erosion then hillslope diffusion\n", - " spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, dt)\n", - " diff_erosion = diffusion_eroder.erode(uplifted_elevation - spl_erosion, dt)\n", - " \n", - " # update topography\n", - " elevation = uplifted_elevation - spl_erosion - diff_erosion\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Plot Outputs and Other Diagnostics\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Topographic elevation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(8, 8))\n", - "plt.imshow(elevation)\n", - "plt.colorbar();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Drainage area (log)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(8, 8))\n", - "plt.imshow(np.log(drainage_area), cmap=plt.cm.Blues)\n", - "plt.colorbar();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Drainage basins" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "colors = [(1,1,1)] + [(random(),random(),random()) for i in range(255)]\n", - "rnd_cm = matplotlib.colors.LinearSegmentedColormap.from_list('new_map', colors, N=256)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(7, 7))\n", - "plt.imshow(flow_graph.basins(), cmap=rnd_cm);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "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.9.10" - } - }, - "nbformat": 4, - "nbformat_minor": 4 + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Inner Base Levels [Py]\n", + "\n", + "An example very similar to {doc}`mountain_py` but setting the base level nodes inside the raster grid instead of on the boundaries. The modeled domain is thus infinite (fully periodic), which is usefull for global (planetary) scale simulations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from random import random, uniform\n", + "\n", + "import fastscapelib as fs\n", + "import numpy as np\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { "tags": ["remove-cell"] }, + "outputs": [], + "source": [ + "# Theme that looks reasonably fine on both dark/light modes\n", + "matplotlib.style.use('Solarize_Light2')\n", + "matplotlib.rcParams['axes.grid'] = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup the Grid, Flow Graph and Eroders\n", + "\n", + "Create a {py:class}`~fastscapelib.RasterGrid` of 201x201 nodes with a total length of 50 km in both y (rows) and x (columns).\n", + "\n", + "Set looped (reflective) boundary conditions at all border nodes. Also set fixed-value \"boundary\" for a given number of base level nodes randomly selected inside the domain." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n_base_levels = 20\n", + "base_level_row = np.random.uniform(1, 200, n_base_levels).astype(\"int\")\n", + "base_level_col = np.random.uniform(1, 200, n_base_levels).astype(\"int\")\n", + "\n", + "base_levels = {\n", + " (i, j): fs.NodeStatus.FIXED_VALUE\n", + " for i, j in zip(base_level_row, base_level_col)\n", + "}\n", + "\n", + "bs = fs.NodeStatus.LOOPED\n", + "grid = fs.RasterGrid.from_length([201, 201], [5e4, 5e4], bs, base_levels)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a {py:class}`~fastscapelib.FlowGraph` object with single direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n", + "\n", + "By default, base level nodes are set from fixed value boundary conditions (random inner nodes in this example)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": + ["flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Setup eroder classes (bedrock channel + hillslope) with a given set of parameter values."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spl_eroder = fs.SPLEroder(\n", + " flow_graph,\n", + " k_coef=2e-4,\n", + " area_exp=0.4,\n", + " slope_exp=1,\n", + " tolerance=1e-5,\n", + ")\n", + "\n", + "diffusion_eroder = fs.DiffusionADIEroder(grid, 0.01)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Setup Initial Conditions and External Forcing"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Create a flat (+ random perturbations) surface topography as initial conditions. Also initialize the array for drainage area."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rng = np.random.Generator(np.random.PCG64(1234))\n", + "\n", + "init_elevation = rng.uniform(0, 1, size=grid.shape)\n", + "\n", + "elevation = init_elevation\n", + "drainage_area = np.empty_like(elevation)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Set upflit rate as uniform (fixed value) within the domain and to zero at all grid boundaries."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "uplift_rate = np.where(\n", + " grid.nodes_status() == fs.NodeStatus.FIXED_VALUE.value, 0, 1e-3\n", + ")" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["## Run the Model"] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Run the model for a few dozens of time steps (total simulation time: 1M years)."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = 2e4\n", + "nsteps = 50\n", + "\n", + "for step in range(nsteps):\n", + " # uplift (no uplift at fixed elevation boundaries)\n", + " uplifted_elevation = elevation + dt * uplift_rate\n", + " \n", + " # flow routing\n", + " filled_elevation = flow_graph.update_routes(uplifted_elevation)\n", + " \n", + " # flow accumulation (drainage area)\n", + " flow_graph.accumulate(drainage_area, 1.0)\n", + " \n", + " # apply channel erosion then hillslope diffusion\n", + " spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, dt)\n", + " diff_erosion = diffusion_eroder.erode(uplifted_elevation - spl_erosion, dt)\n", + " \n", + " # update topography\n", + " elevation = uplifted_elevation - spl_erosion - diff_erosion\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Plot Outputs and Other Diagnostics\n"] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Topographic elevation"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(8, 8))\n", + "plt.imshow(elevation)\n", + "plt.colorbar();" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Drainage area (log)"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(8, 8))\n", + "plt.imshow(np.log(drainage_area), cmap=plt.cm.Blues)\n", + "plt.colorbar();" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Drainage basins"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "colors = [(1,1,1)] + [(random(),random(),random()) for i in range(255)]\n", + "rnd_cm = matplotlib.colors.LinearSegmentedColormap.from_list('new_map', colors, N=256)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(7, 7))\n", + "plt.imshow(flow_graph.basins(), cmap=rnd_cm);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.9.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 } diff --git a/doc/source/examples/mountain_py.ipynb b/doc/source/examples/mountain_py.ipynb index cdf4d4af..cfce7a75 100644 --- a/doc/source/examples/mountain_py.ipynb +++ b/doc/source/examples/mountain_py.ipynb @@ -1,300 +1,256 @@ { "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Mountain [Py]\n", - "\n", - "A simple example simulating the formation of an orogenic belt (on a 2D raster grid) from a nearly-flat surface under the action of block-uplift vs. bedrock channel and hillslope erosion. See also {doc}`mountain_cpp`.\n", - "\n", - "The local rate of elevation change, {math}`\\partial h/\\partial t`, is determined by the balance between uplift (uniform in space and time) {math}`U` and erosion (stream-power law for river erosion and linear diffusion for hillslope erosion).\n", - "\n", - "```{math}\n", - ":label: eq_mountain\n", - "\\frac{\\partial h}{\\partial t} = U - K_f A^m (\\nabla h)^n + K_D \\nabla^2 h\n", - "```\n", - "\n", - "where {math}`A` is the drainage area (i.e., the total upslope area from which the water flow converges) and {math}`K_f`, {math}`m`, {math}`n` and {math}`K_D` are parameters (uniform in space and time).\n", - "\n", - "The initial topography is nearly flat with small random perturbations and elevation remains fixed at the boundaries of the spatial domain." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from random import random\n", - "\n", - "import fastscapelib as fs\n", - "import numpy as np\n", - "import matplotlib\n", - "import matplotlib.pyplot as plt" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "remove-cell" - ] - }, - "outputs": [], - "source": [ - "# Theme that looks reasonably fine on both dark/light modes\n", - "matplotlib.style.use('Solarize_Light2')\n", - "matplotlib.rcParams['axes.grid'] = False" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup the Grid, Flow Graph and Eroders\n", - "\n", - "Create a {py:class}`~fastscapelib.RasterGrid` of 201x301 nodes with a total length of 50 km in y (rows) and 75 km in x (columns).\n", - "\n", - "Set fixed value boundary conditions at all border nodes." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "grid = fs.RasterGrid.from_length([201, 301], [5e4, 7.5e4], fs.NodeStatus.FIXED_VALUE)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create a {py:class}`~fastscapelib.FlowGraph` object with single direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n", - "\n", - "By default, base level nodes are set from fixed value boundary conditions (all border nodes in this example)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setup eroder classes (bedrock channel + hillslope) with a given set of parameter values." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "spl_eroder = fs.SPLEroder(\n", - " flow_graph,\n", - " k_coef=2e-4,\n", - " area_exp=0.4,\n", - " slope_exp=1,\n", - " tolerance=1e-5,\n", - ")\n", - "\n", - "diffusion_eroder = fs.DiffusionADIEroder(grid, 0.01)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup Initial Conditions and External Forcing" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create a flat (+ random perturbations) surface topography as initial conditions. Also initialize the array for drainage area." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "rng = np.random.Generator(np.random.PCG64(1234))\n", - "\n", - "init_elevation = rng.uniform(0, 1, size=grid.shape)\n", - "\n", - "elevation = init_elevation\n", - "drainage_area = np.empty_like(elevation)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Set upflit rate as uniform (fixed value) within the domain and to zero at all grid boundaries." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "uplift_rate = np.full_like(elevation, 1e-3)\n", - "uplift_rate[[0, -1], :] = 0.\n", - "uplift_rate[:, [0, -1]] = 0." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run the Model" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Run the model for a few dozens of time steps (total simulation time: 1M years)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dt = 2e4\n", - "nsteps = 50\n", - "\n", - "for step in range(nsteps):\n", - " # uplift (no uplift at fixed elevation boundaries)\n", - " uplifted_elevation = elevation + dt * uplift_rate\n", - " \n", - " # flow routing\n", - " filled_elevation = flow_graph.update_routes(uplifted_elevation)\n", - " \n", - " # flow accumulation (drainage area)\n", - " flow_graph.accumulate(drainage_area, 1.0)\n", - " \n", - " # apply channel erosion then hillslope diffusion\n", - " spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, dt)\n", - " diff_erosion = diffusion_eroder.erode(uplifted_elevation - spl_erosion, dt)\n", - " \n", - " # update topography\n", - " elevation = uplifted_elevation - spl_erosion - diff_erosion\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Plot Outputs and Other Diagnostics\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Topographic elevation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(12, 6))\n", - "plt.imshow(elevation)\n", - "plt.colorbar();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Drainage area (log)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(12, 6))\n", - "plt.imshow(np.log(drainage_area), cmap=plt.cm.Blues)\n", - "plt.colorbar();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "- Drainage basins" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "colors = [(1,1,1)] + [(random(),random(),random()) for i in range(255)]\n", - "rnd_cm = matplotlib.colors.LinearSegmentedColormap.from_list('new_map', colors, N=256)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(12, 6))\n", - "plt.imshow(flow_graph.basins(), cmap=rnd_cm);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "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.3" - } - }, - "nbformat": 4, - "nbformat_minor": 4 + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Mountain [Py]\n", + "\n", + "A simple example simulating the formation of an orogenic belt (on a 2D raster grid) from a nearly-flat surface under the action of block-uplift vs. bedrock channel and hillslope erosion. See also {doc}`mountain_cpp`.\n", + "\n", + "The local rate of elevation change, {math}`\\partial h/\\partial t`, is determined by the balance between uplift (uniform in space and time) {math}`U` and erosion (stream-power law for river erosion and linear diffusion for hillslope erosion).\n", + "\n", + "```{math}\n", + ":label: eq_mountain\n", + "\\frac{\\partial h}{\\partial t} = U - K_f A^m (\\nabla h)^n + K_D \\nabla^2 h\n", + "```\n", + "\n", + "where {math}`A` is the drainage area (i.e., the total upslope area from which the water flow converges) and {math}`K_f`, {math}`m`, {math}`n` and {math}`K_D` are parameters (uniform in space and time).\n", + "\n", + "The initial topography is nearly flat with small random perturbations and elevation remains fixed at the boundaries of the spatial domain." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from random import random\n", + "\n", + "import fastscapelib as fs\n", + "import numpy as np\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { "tags": ["remove-cell"] }, + "outputs": [], + "source": [ + "# Theme that looks reasonably fine on both dark/light modes\n", + "matplotlib.style.use('Solarize_Light2')\n", + "matplotlib.rcParams['axes.grid'] = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup the Grid, Flow Graph and Eroders\n", + "\n", + "Create a {py:class}`~fastscapelib.RasterGrid` of 201x301 nodes with a total length of 50 km in y (rows) and 75 km in x (columns).\n", + "\n", + "Set fixed value boundary conditions at all border nodes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": + ["grid = fs.RasterGrid.from_length([201, 301], [5e4, 7.5e4], fs.NodeStatus.FIXED_VALUE)"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a {py:class}`~fastscapelib.FlowGraph` object with single direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n", + "\n", + "By default, base level nodes are set from fixed value boundary conditions (all border nodes in this example)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": + ["flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Setup eroder classes (bedrock channel + hillslope) with a given set of parameter values."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spl_eroder = fs.SPLEroder(\n", + " flow_graph,\n", + " k_coef=2e-4,\n", + " area_exp=0.4,\n", + " slope_exp=1,\n", + " tolerance=1e-5,\n", + ")\n", + "\n", + "diffusion_eroder = fs.DiffusionADIEroder(grid, 0.01)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Setup Initial Conditions and External Forcing"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Create a flat (+ random perturbations) surface topography as initial conditions. Also initialize the array for drainage area."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rng = np.random.Generator(np.random.PCG64(1234))\n", + "\n", + "init_elevation = rng.uniform(0, 1, size=grid.shape)\n", + "\n", + "elevation = init_elevation\n", + "drainage_area = np.empty_like(elevation)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Set upflit rate as uniform (fixed value) within the domain and to zero at all grid boundaries."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "uplift_rate = np.full_like(elevation, 1e-3)\n", + "uplift_rate[[0, -1], :] = 0.\n", + "uplift_rate[:, [0, -1]] = 0." + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["## Run the Model"] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Run the model for a few dozens of time steps (total simulation time: 1M years)."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = 2e4\n", + "nsteps = 50\n", + "\n", + "for step in range(nsteps):\n", + " # uplift (no uplift at fixed elevation boundaries)\n", + " uplifted_elevation = elevation + dt * uplift_rate\n", + " \n", + " # flow routing\n", + " filled_elevation = flow_graph.update_routes(uplifted_elevation)\n", + " \n", + " # flow accumulation (drainage area)\n", + " flow_graph.accumulate(drainage_area, 1.0)\n", + " \n", + " # apply channel erosion then hillslope diffusion\n", + " spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, dt)\n", + " diff_erosion = diffusion_eroder.erode(uplifted_elevation - spl_erosion, dt)\n", + " \n", + " # update topography\n", + " elevation = uplifted_elevation - spl_erosion - diff_erosion\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Plot Outputs and Other Diagnostics\n"] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Topographic elevation"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 6))\n", + "plt.imshow(elevation)\n", + "plt.colorbar();" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Drainage area (log)"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 6))\n", + "plt.imshow(np.log(drainage_area), cmap=plt.cm.Blues)\n", + "plt.colorbar();" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["- Drainage basins"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "colors = [(1,1,1)] + [(random(),random(),random()) for i in range(255)]\n", + "rnd_cm = matplotlib.colors.LinearSegmentedColormap.from_list('new_map', colors, N=256)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 6))\n", + "plt.imshow(flow_graph.basins(), cmap=rnd_cm);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 } diff --git a/doc/source/examples/river_profile_py.ipynb b/doc/source/examples/river_profile_py.ipynb index 0c722f9f..e747f5d1 100644 --- a/doc/source/examples/river_profile_py.ipynb +++ b/doc/source/examples/river_profile_py.ipynb @@ -1,271 +1,230 @@ { "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# River Profile [Py]\n", - "\n", - "A simple example simulating the evolution of a river longitudinal profile under the action of block-uplift vs. bedrock channel erosion, starting from a very gentle channel slope. See also {doc}`river_profile_cpp`.\n", - "\n", - "The local rate of elevation change, {math}`\\partial h/\\partial t`, is determined by the balance between uplift (uniform in space and time) {math}`U` and bedrock channel erosion modeled using the Stream-Power Law (SPL).\n", - "\n", - "```{math}\n", - ":label: eq_river\n", - "\\frac{\\partial h}{\\partial t} = U - K_f A^m \\left(\\frac{\\partial h}{\\partial x}\\right)^n\n", - "```\n", - "\n", - "where {math}`A` is the drainage area and {math}`K_f`, {math}`m`, {math}`n` are SPL parameters (uniform in space and time)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import fastscapelib as fs\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "remove-cell" - ] - }, - "outputs": [], - "source": [ - "# Theme that looks reasonably fine on both dark/light modes\n", - "import matplotlib\n", - "matplotlib.style.use('Solarize_Light2')\n", - "matplotlib.rcParams['axes.grid'] = False" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup the Grid, Flow Graph and Eroders\n", - "\n", - "Create a {py:class}`~fastscapelib.ProfileGrid` of 101 nodes each uniformly spaced by 300 meters.\n", - "\n", - "Set fixed-value boundary conditions on the left side and free (core) conditions on the right side of the profile." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "grid = fs.ProfileGrid(101, 300.0, [fs.NodeStatus.FIXED_VALUE, fs.NodeStatus.CORE])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Set x-coordinate values along the profile grid, which correspond to the distance from the ridge top. Start at `x0 > 0` so that the channel head effectively starts at some distance away from the ridge top." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x0 = 300.0\n", - "length = (grid.size - 1) * grid.spacing\n", - "x = np.linspace(x0 + length, x0, num=grid.size)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create a {py:class}`~fastscapelib.FlowGraph` object with single direction flow routing (trivial case for a 1-dimensional profile).\n", - "\n", - "By default, base level nodes are set from fixed value boundary conditions (left node in this example)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter()])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Compute drainage area along the profile using [Hack's Law](https://en.wikipedia.org/wiki/Hack%27s_law)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "hack_coef = 6.69\n", - "hack_exp = 1.67\n", - "\n", - "drainage_area = hack_coef * x**hack_exp" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setup the SPL eroder class with a given set of parameter values." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "spl_eroder = fs.SPLEroder(\n", - " flow_graph,\n", - " k_coef=1e-4,\n", - " area_exp=0.5,\n", - " slope_exp=1,\n", - " tolerance=1e-5,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup Initial Conditions and External Forcing" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Create an initial river profile with a very gentle, uniform slope." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "init_elevation = (length + x0 - x) * 1e-4\n", - "elevation = init_elevation" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Set upflit rate as uniform (fixed value) along the profile except at the base level node (zero)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "uplift_rate = np.full_like(x, 1e-3)\n", - "uplift_rate[0] = 0." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run the Model" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Run the model for a few thousands of time steps. At the middle of the simulation, change the value of the SPL coefficient to simulate a change in erodibility (e.g., climate change). This should create a knicpoint in the resulting river profile." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dt = 1e2\n", - "nsteps = 4000\n", - "\n", - "# flow routing is required for solving SPL, although in this example\n", - "# it can be computed out of the time step loop since the\n", - "# flow path (single river) will not change during the simulation\n", - "# (not the case, e.g., if there were closed depressions!)\n", - "flow_graph.update_routes(elevation)\n", - "\n", - "\n", - "for step in range(nsteps):\n", - " # uplift (no uplift at fixed elevation boundaries)\n", - " uplifted_elevation = elevation + dt * uplift_rate\n", - " \n", - " # abrubpt change in erodibility\n", - " if step == 3000:\n", - " spl_eroder.k_coef /= 4\n", - "\n", - " # apply channel erosion\n", - " spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, dt)\n", - " \n", - " # update topography\n", - " elevation = uplifted_elevation - spl_erosion\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Plot the River Profile" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(12, 6))\n", - "ax.plot(x, elevation)\n", - "plt.setp(ax, xlabel=\"x\", ylabel=\"elevation\", xlim=(x[0] + 300.0, 0));" - ] - } - ], - "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.9.10" - } - }, - "nbformat": 4, - "nbformat_minor": 4 + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# River Profile [Py]\n", + "\n", + "A simple example simulating the evolution of a river longitudinal profile under the action of block-uplift vs. bedrock channel erosion, starting from a very gentle channel slope. See also {doc}`river_profile_cpp`.\n", + "\n", + "The local rate of elevation change, {math}`\\partial h/\\partial t`, is determined by the balance between uplift (uniform in space and time) {math}`U` and bedrock channel erosion modeled using the Stream-Power Law (SPL).\n", + "\n", + "```{math}\n", + ":label: eq_river\n", + "\\frac{\\partial h}{\\partial t} = U - K_f A^m \\left(\\frac{\\partial h}{\\partial x}\\right)^n\n", + "```\n", + "\n", + "where {math}`A` is the drainage area and {math}`K_f`, {math}`m`, {math}`n` are SPL parameters (uniform in space and time)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import fastscapelib as fs\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { "tags": ["remove-cell"] }, + "outputs": [], + "source": [ + "# Theme that looks reasonably fine on both dark/light modes\n", + "import matplotlib\n", + "matplotlib.style.use('Solarize_Light2')\n", + "matplotlib.rcParams['axes.grid'] = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup the Grid, Flow Graph and Eroders\n", + "\n", + "Create a {py:class}`~fastscapelib.ProfileGrid` of 101 nodes each uniformly spaced by 300 meters.\n", + "\n", + "Set fixed-value boundary conditions on the left side and free (core) conditions on the right side of the profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": + ["grid = fs.ProfileGrid(101, 300.0, [fs.NodeStatus.FIXED_VALUE, fs.NodeStatus.CORE])"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Set x-coordinate values along the profile grid, which correspond to the distance from the ridge top. Start at `x0 > 0` so that the channel head effectively starts at some distance away from the ridge top."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x0 = 300.0\n", + "length = (grid.size - 1) * grid.spacing\n", + "x = np.linspace(x0 + length, x0, num=grid.size)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a {py:class}`~fastscapelib.FlowGraph` object with single direction flow routing (trivial case for a 1-dimensional profile).\n", + "\n", + "By default, base level nodes are set from fixed value boundary conditions (left node in this example)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": ["flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter()])"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Compute drainage area along the profile using [Hack's Law](https://en.wikipedia.org/wiki/Hack%27s_law)."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "hack_coef = 6.69\n", + "hack_exp = 1.67\n", + "\n", + "drainage_area = hack_coef * x**hack_exp" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["Setup the SPL eroder class with a given set of parameter values."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spl_eroder = fs.SPLEroder(\n", + " flow_graph,\n", + " k_coef=1e-4,\n", + " area_exp=0.5,\n", + " slope_exp=1,\n", + " tolerance=1e-5,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["## Setup Initial Conditions and External Forcing"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": ["Create an initial river profile with a very gentle, uniform slope."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": ["init_elevation = (length + x0 - x) * 1e-4\n", "elevation = init_elevation"] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Set upflit rate as uniform (fixed value) along the profile except at the base level node (zero)."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": ["uplift_rate = np.full_like(x, 1e-3)\n", "uplift_rate[0] = 0."] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["## Run the Model"] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": + ["Run the model for a few thousands of time steps. At the middle of the simulation, change the value of the SPL coefficient to simulate a change in erodibility (e.g., climate change). This should create a knicpoint in the resulting river profile."] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dt = 1e2\n", + "nsteps = 4000\n", + "\n", + "# flow routing is required for solving SPL, although in this example\n", + "# it can be computed out of the time step loop since the\n", + "# flow path (single river) will not change during the simulation\n", + "# (not the case, e.g., if there were closed depressions!)\n", + "flow_graph.update_routes(elevation)\n", + "\n", + "\n", + "for step in range(nsteps):\n", + " # uplift (no uplift at fixed elevation boundaries)\n", + " uplifted_elevation = elevation + dt * uplift_rate\n", + " \n", + " # abrubpt change in erodibility\n", + " if step == 3000:\n", + " spl_eroder.k_coef /= 4\n", + "\n", + " # apply channel erosion\n", + " spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, dt)\n", + " \n", + " # update topography\n", + " elevation = uplifted_elevation - spl_erosion\n" + ] + }, + { "cell_type": "markdown", "metadata": {}, "source": ["## Plot the River Profile"] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 6))\n", + "ax.plot(x, elevation)\n", + "plt.setp(ax, xlabel=\"x\", ylabel=\"elevation\", xlim=(x[0] + 300.0, 0));" + ] + } + ], + "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.9.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 } diff --git a/examples/eroder_kernel.ipynb b/examples/eroder_kernel.ipynb index 6107e0d2..580ba536 100644 --- a/examples/eroder_kernel.ipynb +++ b/examples/eroder_kernel.ipynb @@ -1,275 +1,259 @@ { "cells": [ - { - "cell_type": "markdown", - "id": "076b1ff5-7333-438b-a25f-735434af0543", - "metadata": {}, - "source": [ - "# Fastscapelib `NumbaEroderType` example" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1c79bc75-4c5b-471f-aa79-dd44e1069f1b", - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numba as nb\n", - "import numpy as np\n", - "\n", - "from fastscapelib.eroders import FlowKernelEroder\n", - "import fastscapelib as fs\n", - "\n", - "\n", - "grid = fs.RasterGrid.from_length([201, 301], [5e4, 7.5e4], fs.NodeStatus.FIXED_VALUE)\n", - "# flow_graph = fs.FlowGraph(grid, [fs.MultiFlowRouter()])\n", - "flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])\n", - "\n", - "rng = np.random.Generator(np.random.PCG64(1234))\n", - "\n", - "init_elevation = rng.uniform(0, 5, size=grid.shape)\n", - "drainage_area = np.empty_like(init_elevation)\n", - "\n", - "uplift_rate = np.full_like(init_elevation, 1e-3)\n", - "uplift_rate[[0, -1], :] = 0.0\n", - "uplift_rate[:, [0, -1]] = 0.0\n", - "flow_graph.update_routes(init_elevation)\n", - "\n", - "class NumbaSplEroder(FlowKernelEroder):\n", - " spec = dict(\n", - " elevation=nb.float64[::1],\n", - " erosion=nb.float64[::1],\n", - " drainage_area=nb.float64[::1],\n", - " k_coef=(nb.float64, 2e-4),\n", - " area_exp=(nb.float64,0.4),\n", - " slope_exp=(nb.float64, 1.),\n", - " dt=(nb.float64, 2e4),\n", - " )\n", - "\n", - " outputs = [\"erosion\"]\n", - " \n", - " apply_dir=fs.FlowGraphTraversalDir.BREADTH_UPSTREAM\n", - " \n", - " def __init__(self, *args, **kwargs):\n", - " super().__init__(*args, **kwargs)\n", - " self.kernel_data.bind(erosion=np.zeros(flow_graph.size))\n", - " \n", - " def erode(self, elevation, drainage_area, dt):\n", - " self.set_kernel_data(elevation=elevation, drainage_area=drainage_area, dt=dt)\n", - " self.kernel_data.erosion.fill(0.)\n", - " super().erode()\n", - " return self.kernel_data.erosion\n", - " \n", - " @staticmethod\n", - " def eroder_kernel(node):\n", - " dt = node.dt\n", - "\n", - " r_count = node.receivers.count\n", - " if r_count == 1 and node.receivers.distance[0] == 0.0:\n", - " return\n", - " \n", - " elevation_flooded = np.finfo(np.double).max\n", - " \n", - " for r in range(r_count):\n", - " irec_elevation_next = node.receivers.elevation[r] - node.receivers.erosion[r]\n", - " \n", - " if irec_elevation_next < elevation_flooded:\n", - " elevation_flooded = irec_elevation_next\n", - " \n", - " if node.elevation <= elevation_flooded:\n", - " return\n", - " \n", - " eq_num = node.elevation\n", - " eq_den = 1.0\n", - " \n", - " for r in range(r_count):\n", - " irec_elevation = node.receivers.elevation[r]\n", - " irec_elevation_next = irec_elevation - node.receivers.erosion[r]\n", - " \n", - " if irec_elevation > node.elevation:\n", - " continue\n", - " \n", - " irec_weight = node.receivers.weight[r]\n", - " irec_distance = node.receivers.distance[r]\n", - " \n", - " factor = (\n", - " node.k_coef * node.dt * np.power(node.drainage_area * irec_weight, node.area_exp)\n", - " )\n", - " factor /= irec_distance\n", - " eq_num += factor * irec_elevation_next\n", - " eq_den += factor\n", - " \n", - " elevation_updated = eq_num / eq_den\n", - " \n", - " if elevation_updated < elevation_flooded:\n", - " elevation_updated = elevation_flooded + np.finfo(np.double).tiny\n", - " \n", - " node.erosion = node.elevation - elevation_updated\n", - " \n", - "elevation = init_elevation.ravel().copy()\n", - "erosion = np.zeros(flow_graph.size)\n", - "drainage_area = np.ones(flow_graph.size)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "17991b9f-69fb-44e3-abed-c46af557a8d1", - "metadata": {}, - "outputs": [], - "source": [ - "dt = 2e4\n", - "nsteps = 50\n", - "uplift = dt * uplift_rate.ravel()\n", - "\n", - "spl_eroder = NumbaSplEroder(flow_graph, max_receivers=10)\n", - "\n", - "fs_spl_eroder = fs.SPLEroder(\n", - " flow_graph,\n", - " k_coef=2e-4,\n", - " area_exp=0.4,\n", - " slope_exp=1.,\n", - " tolerance=1e-5,\n", - ")\n", - "\n", - "def run_simulation(eroder):\n", - " elevation = init_elevation.copy().ravel()\n", - " drainage_area = np.empty_like(elevation)\n", - "\n", - " for step in range(nsteps):\n", - " # uplift (no uplift at fixed elevation boundaries)\n", - " uplifted_elevation = elevation + uplift\n", - " \n", - " # flow routing\n", - " flow_graph.update_routes(uplifted_elevation)\n", - "\n", - " # flow accumulation (drainage area)\n", - " flow_graph.accumulate(drainage_area, 1.0)\n", - "\n", - " erosion = eroder.erode(uplifted_elevation, drainage_area, dt)\n", - "\n", - " # update topography\n", - " elevation = uplifted_elevation - erosion.ravel()\n", - "\n", - " return elevation.reshape(grid.shape)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a1186cf5-29c8-4107-82eb-eb0aebdab31c", - "metadata": {}, - "outputs": [], - "source": [ - "spl_eroder.erode(elevation=elevation, drainage_area=drainage_area, dt=dt);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c80f6143-445e-4d8f-8520-c54b7b86e7bb", - "metadata": {}, - "outputs": [], - "source": [ - "for i in [1, 2, 4, 8]:\n", - " spl_eroder.n_threads = i\n", - " %timeit -r 10 -n 10 spl_eroder.erode(elevation=elevation, drainage_area=drainage_area, dt=dt)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "90de6ce1-4ba7-4225-a60c-6f957b2c2691", - "metadata": {}, - "outputs": [], - "source": [ - "for i in [1, 2, 4, 8, 16]:\n", - " spl_eroder.n_threads = i\n", - " %timeit -r 1 -n 1 run_simulation(spl_eroder)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b25dde09-f598-4a85-8baa-796556d31b7d", - "metadata": {}, - "outputs": [], - "source": [ - "# efficiency of the parallel execution depends to the kernel size, the current flow graph and the kernel execution order\n", - "for i in [1, 2, 4, 8]:\n", - " spl_eroder.n_threads = i\n", - " %timeit -r 10 -n 10 spl_eroder.erode(elevation=elevation, drainage_area=drainage_area, dt=dt)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f7089e1f-ca26-42a9-89ee-e2711d5216d9", - "metadata": {}, - "outputs": [], - "source": [ - "%timeit -r 1 -n 1 run_simulation(fs_spl_eroder)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8db53d5a-a17d-4079-86d0-bed75bc1d999", - "metadata": {}, - "outputs": [], - "source": [ - "%timeit -r 10 -n 10 fs_spl_eroder.erode(elevation, drainage_area, dt)" - ] - }, - { - "cell_type": "markdown", - "id": "a50faaa5-8c71-4d7c-acce-faa733f73c49", - "metadata": {}, - "source": [ - "### Benchmark Results" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7ab15183-51e6-4072-8e69-a06c24d6757c", - "metadata": {}, - "outputs": [], - "source": [ - "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(20, 10))\n", - "\n", - "spl_eroder._kernel.kernel.n_threads = 1\n", - "axes[0, 0].imshow(run_simulation(spl_eroder).reshape(grid.shape))\n", - "spl_eroder._kernel.kernel.n_threads = 4\n", - "axes[0, 1].imshow(run_simulation(spl_eroder).reshape(grid.shape))\n", - "spl_eroder._kernel.kernel.n_threads = 8\n", - "axes[1, 0].imshow(run_simulation(spl_eroder).reshape(grid.shape))\n", - "\n", - "axes[1, 1].imshow(run_simulation(fs_spl_eroder));" - ] - } - ], - "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.3" - } - }, - "nbformat": 4, - "nbformat_minor": 5 + { + "cell_type": "markdown", + "id": "076b1ff5-7333-438b-a25f-735434af0543", + "metadata": {}, + "source": ["# Fastscapelib `NumbaEroderType` example"] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c79bc75-4c5b-471f-aa79-dd44e1069f1b", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numba as nb\n", + "import numpy as np\n", + "\n", + "from fastscapelib.eroders import FlowKernelEroder\n", + "import fastscapelib as fs\n", + "\n", + "\n", + "grid = fs.RasterGrid.from_length([201, 301], [5e4, 7.5e4], fs.NodeStatus.FIXED_VALUE)\n", + "# flow_graph = fs.FlowGraph(grid, [fs.MultiFlowRouter()])\n", + "flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])\n", + "\n", + "rng = np.random.Generator(np.random.PCG64(1234))\n", + "\n", + "init_elevation = rng.uniform(0, 5, size=grid.shape)\n", + "drainage_area = np.empty_like(init_elevation)\n", + "\n", + "uplift_rate = np.full_like(init_elevation, 1e-3)\n", + "uplift_rate[[0, -1], :] = 0.0\n", + "uplift_rate[:, [0, -1]] = 0.0\n", + "flow_graph.update_routes(init_elevation)\n", + "\n", + "class NumbaSplEroder(FlowKernelEroder):\n", + " spec = dict(\n", + " elevation=nb.float64[::1],\n", + " erosion=nb.float64[::1],\n", + " drainage_area=nb.float64[::1],\n", + " k_coef=(nb.float64, 2e-4),\n", + " area_exp=(nb.float64,0.4),\n", + " slope_exp=(nb.float64, 1.),\n", + " dt=(nb.float64, 2e4),\n", + " )\n", + "\n", + " outputs = [\"erosion\"]\n", + " \n", + " apply_dir=fs.FlowGraphTraversalDir.BREADTH_UPSTREAM\n", + " \n", + " def __init__(self, *args, **kwargs):\n", + " super().__init__(*args, **kwargs)\n", + " self.kernel_data.bind(erosion=np.zeros(flow_graph.size))\n", + " \n", + " def erode(self, elevation, drainage_area, dt):\n", + " self.set_kernel_data(elevation=elevation, drainage_area=drainage_area, dt=dt)\n", + " self.kernel_data.erosion.fill(0.)\n", + " super().erode()\n", + " return self.kernel_data.erosion\n", + " \n", + " @staticmethod\n", + " def eroder_kernel(node):\n", + " dt = node.dt\n", + "\n", + " r_count = node.receivers.count\n", + " if r_count == 1 and node.receivers.distance[0] == 0.0:\n", + " return\n", + " \n", + " elevation_flooded = np.finfo(np.double).max\n", + " \n", + " for r in range(r_count):\n", + " irec_elevation_next = node.receivers.elevation[r] - node.receivers.erosion[r]\n", + " \n", + " if irec_elevation_next < elevation_flooded:\n", + " elevation_flooded = irec_elevation_next\n", + " \n", + " if node.elevation <= elevation_flooded:\n", + " return\n", + " \n", + " eq_num = node.elevation\n", + " eq_den = 1.0\n", + " \n", + " for r in range(r_count):\n", + " irec_elevation = node.receivers.elevation[r]\n", + " irec_elevation_next = irec_elevation - node.receivers.erosion[r]\n", + " \n", + " if irec_elevation > node.elevation:\n", + " continue\n", + " \n", + " irec_weight = node.receivers.weight[r]\n", + " irec_distance = node.receivers.distance[r]\n", + " \n", + " factor = (\n", + " node.k_coef * node.dt * np.power(node.drainage_area * irec_weight, node.area_exp)\n", + " )\n", + " factor /= irec_distance\n", + " eq_num += factor * irec_elevation_next\n", + " eq_den += factor\n", + " \n", + " elevation_updated = eq_num / eq_den\n", + " \n", + " if elevation_updated < elevation_flooded:\n", + " elevation_updated = elevation_flooded + np.finfo(np.double).tiny\n", + " \n", + " node.erosion = node.elevation - elevation_updated\n", + " \n", + "elevation = init_elevation.ravel().copy()\n", + "erosion = np.zeros(flow_graph.size)\n", + "drainage_area = np.ones(flow_graph.size)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17991b9f-69fb-44e3-abed-c46af557a8d1", + "metadata": {}, + "outputs": [], + "source": [ + "dt = 2e4\n", + "nsteps = 50\n", + "uplift = dt * uplift_rate.ravel()\n", + "\n", + "spl_eroder = NumbaSplEroder(flow_graph, max_receivers=10)\n", + "\n", + "fs_spl_eroder = fs.SPLEroder(\n", + " flow_graph,\n", + " k_coef=2e-4,\n", + " area_exp=0.4,\n", + " slope_exp=1.,\n", + " tolerance=1e-5,\n", + ")\n", + "\n", + "def run_simulation(eroder):\n", + " elevation = init_elevation.copy().ravel()\n", + " drainage_area = np.empty_like(elevation)\n", + "\n", + " for step in range(nsteps):\n", + " # uplift (no uplift at fixed elevation boundaries)\n", + " uplifted_elevation = elevation + uplift\n", + " \n", + " # flow routing\n", + " flow_graph.update_routes(uplifted_elevation)\n", + "\n", + " # flow accumulation (drainage area)\n", + " flow_graph.accumulate(drainage_area, 1.0)\n", + "\n", + " erosion = eroder.erode(uplifted_elevation, drainage_area, dt)\n", + "\n", + " # update topography\n", + " elevation = uplifted_elevation - erosion.ravel()\n", + "\n", + " return elevation.reshape(grid.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1186cf5-29c8-4107-82eb-eb0aebdab31c", + "metadata": {}, + "outputs": [], + "source": ["spl_eroder.erode(elevation=elevation, drainage_area=drainage_area, dt=dt);"] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c80f6143-445e-4d8f-8520-c54b7b86e7bb", + "metadata": {}, + "outputs": [], + "source": [ + "for i in [1, 2, 4, 8]:\n", + " spl_eroder.n_threads = i\n", + " %timeit -r 10 -n 10 spl_eroder.erode(elevation=elevation, drainage_area=drainage_area, dt=dt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90de6ce1-4ba7-4225-a60c-6f957b2c2691", + "metadata": {}, + "outputs": [], + "source": [ + "for i in [1, 2, 4, 8, 16]:\n", + " spl_eroder.n_threads = i\n", + " %timeit -r 1 -n 1 run_simulation(spl_eroder)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b25dde09-f598-4a85-8baa-796556d31b7d", + "metadata": {}, + "outputs": [], + "source": [ + "# efficiency of the parallel execution depends to the kernel size, the current flow graph and the kernel execution order\n", + "for i in [1, 2, 4, 8]:\n", + " spl_eroder.n_threads = i\n", + " %timeit -r 10 -n 10 spl_eroder.erode(elevation=elevation, drainage_area=drainage_area, dt=dt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7089e1f-ca26-42a9-89ee-e2711d5216d9", + "metadata": {}, + "outputs": [], + "source": ["%timeit -r 1 -n 1 run_simulation(fs_spl_eroder)"] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8db53d5a-a17d-4079-86d0-bed75bc1d999", + "metadata": {}, + "outputs": [], + "source": ["%timeit -r 10 -n 10 fs_spl_eroder.erode(elevation, drainage_area, dt)"] + }, + { + "cell_type": "markdown", + "id": "a50faaa5-8c71-4d7c-acce-faa733f73c49", + "metadata": {}, + "source": ["### Benchmark Results"] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ab15183-51e6-4072-8e69-a06c24d6757c", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(20, 10))\n", + "\n", + "spl_eroder._kernel.kernel.n_threads = 1\n", + "axes[0, 0].imshow(run_simulation(spl_eroder).reshape(grid.shape))\n", + "spl_eroder._kernel.kernel.n_threads = 4\n", + "axes[0, 1].imshow(run_simulation(spl_eroder).reshape(grid.shape))\n", + "spl_eroder._kernel.kernel.n_threads = 8\n", + "axes[1, 0].imshow(run_simulation(spl_eroder).reshape(grid.shape))\n", + "\n", + "axes[1, 1].imshow(run_simulation(fs_spl_eroder));" + ] + } + ], + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 } diff --git a/examples/flow_kernel.ipynb b/examples/flow_kernel.ipynb index f416569d..3b8c78de 100644 --- a/examples/flow_kernel.ipynb +++ b/examples/flow_kernel.ipynb @@ -1,335 +1,319 @@ { "cells": [ - { - "cell_type": "markdown", - "id": "076b1ff5-7333-438b-a25f-735434af0543", - "metadata": {}, - "source": [ - "# Fastscapelib \"user flow kernels\" prototype" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1c79bc75-4c5b-471f-aa79-dd44e1069f1b", - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numba as nb\n", - "import numpy as np\n", - "\n", - "import fastscapelib as fs\n", - "\n", - "grid = fs.RasterGrid.from_length([201, 301], [5e4, 7.5e4], fs.NodeStatus.FIXED_VALUE)\n", - "flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])\n", - "\n", - "rng = np.random.Generator(np.random.PCG64(1234))\n", - "\n", - "init_elevation = rng.uniform(0, 5, size=grid.shape)\n", - "drainage_area = np.empty_like(init_elevation)\n", - "\n", - "uplift_rate = np.full_like(init_elevation, 1e-3)\n", - "uplift_rate[[0, -1], :] = 0.0\n", - "uplift_rate[:, [0, -1]] = 0.0\n", - "flow_graph.update_routes(init_elevation)\n", - "\n", - "\n", - "def kernel_func(node): \n", - " dt = node.dt\n", - " r_count = node.receivers.count\n", - " if r_count == 1 and node.receivers.distance[0] == 0.0:\n", - " return\n", - "\n", - " elevation_flooded = np.finfo(np.double).max\n", - "\n", - " for r in range(r_count):\n", - " irec_elevation_next = node.receivers.elevation[r] - node.receivers.erosion[r]\n", - "\n", - " if irec_elevation_next < elevation_flooded:\n", - " elevation_flooded = irec_elevation_next\n", - "\n", - " if node.elevation <= elevation_flooded:\n", - " return\n", - "\n", - " eq_num = node.elevation\n", - " eq_den = 1.0\n", - "\n", - " for r in range(r_count):\n", - " irec_elevation = node.receivers.elevation[r]\n", - " irec_elevation_next = irec_elevation - node.receivers.erosion[r]\n", - "\n", - " if irec_elevation > node.elevation:\n", - " continue\n", - "\n", - " irec_weight = node.receivers.weight[r]\n", - " irec_distance = node.receivers.distance[r]\n", - "\n", - " factor = (\n", - " node.k_coef * node.dt * np.power(node.drainage_area * irec_weight, node.area_exp)\n", - " )\n", - " factor /= irec_distance\n", - " eq_num += factor * irec_elevation_next\n", - " eq_den += factor\n", - "\n", - " elevation_updated = eq_num / eq_den\n", - "\n", - " if elevation_updated < elevation_flooded:\n", - " elevation_updated = elevation_flooded + np.finfo(np.double).tiny\n", - "\n", - " node.erosion = node.elevation - elevation_updated\n", - "\n", - "elevation = init_elevation.ravel().copy()\n", - "erosion = np.zeros(flow_graph.size)\n", - "drainage_area = np.ones(flow_graph.size)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d4222ec0-82c0-4617-abef-9dfd44e30ad8", - "metadata": {}, - "outputs": [], - "source": [ - "kernel, data = fs.create_flow_kernel(\n", - " flow_graph,\n", - " kernel_func,\n", - " spec=dict(\n", - " elevation=nb.float64[::1],\n", - " erosion=nb.float64[::1],\n", - " drainage_area=nb.float64[::1],\n", - " k_coef=(nb.float64, 2e-4),\n", - " area_exp=(nb.float64, 0.4),\n", - " slope_exp=(nb.float64, 1.),\n", - " dt=(nb.float64, 2e4),\n", - " ),\n", - " outputs=[\"erosion\"],\n", - " max_receivers=1,\n", - " n_threads=1,\n", - " # print_generated_code=True,\n", - " apply_dir=fs.flow.FlowGraphTraversalDir.BREADTH_UPSTREAM\n", - ")\n", - "\n", - "data.bind(\n", - " elevation=elevation,\n", - " erosion=erosion,\n", - " drainage_area=drainage_area,\n", - ")\n", - "\n", - "flow_graph.apply_kernel(kernel, data)\n", - "\n", - "import time\n", - "N_TIMES = 1\n", - "\n", - "t1 = time.time()\n", - "for _ in range(N_TIMES):\n", - " flow_graph.apply_kernel(kernel, data)\n", - "print(round((time.time() - t1) * 1e9 / N_TIMES / grid.size, 1))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "17991b9f-69fb-44e3-abed-c46af557a8d1", - "metadata": {}, - "outputs": [], - "source": [ - "data.dt = dt = 2e4\n", - "nsteps = 50\n", - "uplift = dt * uplift_rate.ravel()\n", - "\n", - "def run_simulation():\n", - " elevation = init_elevation.copy().ravel()\n", - " drainage_area = np.empty_like(elevation)\n", - "\n", - " for step in range(nsteps):\n", - " # uplift (no uplift at fixed elevation boundaries)\n", - " uplifted_elevation = elevation + uplift\n", - " \n", - " # flow routing\n", - " flow_graph.update_routes(uplifted_elevation)\n", - "\n", - " # flow accumulation (drainage area)\n", - " flow_graph.accumulate(drainage_area, 1.0)\n", - "\n", - " data.erosion.fill(0.)\n", - " data.bind(elevation=uplifted_elevation, drainage_area=drainage_area, dt=dt)\n", - " flow_graph.apply_kernel(kernel, data)\n", - "\n", - " # update topography\n", - " elevation = uplifted_elevation - data.erosion\n", - "\n", - " return elevation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c80f6143-445e-4d8f-8520-c54b7b86e7bb", - "metadata": {}, - "outputs": [], - "source": [ - "for i in [1, 2, 4, 8]:\n", - " kernel.kernel.n_threads = i\n", - " %timeit -r 10 -n 10 flow_graph.apply_kernel(kernel, data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "90de6ce1-4ba7-4225-a60c-6f957b2c2691", - "metadata": {}, - "outputs": [], - "source": [ - "for i in [1, 2, 4, 8]:\n", - " kernel.kernel.n_threads = i\n", - " %timeit -r 1 -n 1 run_simulation()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a6598a97-acf5-4a14-a7e6-5310e3688a6d", - "metadata": {}, - "outputs": [], - "source": [ - "# efficiency of the parallel execution depends to the kernel size, the current flow graph and the kernel execution order\n", - "for i in [1, 2, 4, 8]:\n", - " kernel.kernel.n_threads = i\n", - " %timeit -r 10 -n 10 flow_graph.apply_kernel(kernel, data)" - ] - }, - { - "cell_type": "markdown", - "id": "1235bd7e-0dae-465a-ab6f-630b0a6e3cb5", - "metadata": {}, - "source": [ - "### Fastscapelib Python bindings" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "10c5fae2-dbf7-4971-9c16-73429ff9bca0", - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numba as nb\n", - "import numpy as np\n", - "\n", - "import fastscapelib as fs\n", - "\n", - "uplift_rate = np.full_like(init_elevation, 1e-3)\n", - "uplift_rate[[0, -1], :] = 0.\n", - "uplift_rate[:, [0, -1]] = 0.\n", - "\n", - "dt = 2e4\n", - "nsteps = 50\n", - "\n", - "k_coef = 2e-4\n", - "area_exp = 0.4\n", - "slope_exp = 1\n", - "\n", - "flow_graph.update_routes(init_elevation)\n", - "\n", - "spl_eroder = fs.SPLEroder(\n", - " flow_graph,\n", - " k_coef=k_coef,\n", - " area_exp=area_exp,\n", - " slope_exp=slope_exp,\n", - " tolerance=1e-5,\n", - ")\n", - "def run_simulation2():\n", - " elevation = init_elevation.copy()\n", - " drainage_area = np.empty_like(init_elevation)\n", - "\n", - " for step in range(nsteps):\n", - " # uplift (no uplift at fixed elevation boundaries)\n", - " uplifted_elevation = elevation + dt * uplift_rate\n", - " \n", - " # flow routing\n", - " filled_elevation = flow_graph.update_routes(uplifted_elevation)\n", - " \n", - " # flow accumulation (drainage area)\n", - " flow_graph.accumulate(drainage_area, 1.0)\n", - " \n", - " # apply channel erosion (SPL)\n", - " spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, dt)\n", - "\n", - " # update topography\n", - " elevation = uplifted_elevation - spl_erosion\n", - "\n", - " return elevation " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f7089e1f-ca26-42a9-89ee-e2711d5216d9", - "metadata": {}, - "outputs": [], - "source": [ - "%timeit -r 1 -n 1 run_simulation2()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8db53d5a-a17d-4079-86d0-bed75bc1d999", - "metadata": {}, - "outputs": [], - "source": [ - "%timeit -r 10 -n 10 spl_eroder.erode(elevation, drainage_area, dt)" - ] - }, - { - "cell_type": "markdown", - "id": "a50faaa5-8c71-4d7c-acce-faa733f73c49", - "metadata": {}, - "source": [ - "### Benchmark Results" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7ab15183-51e6-4072-8e69-a06c24d6757c", - "metadata": {}, - "outputs": [], - "source": [ - "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(20, 10))\n", - "\n", - "kernel.kernel.n_threads = 1\n", - "axes[0, 0].imshow(run_simulation().reshape(grid.shape))\n", - "kernel.kernel.n_threads = 4\n", - "axes[0, 1].imshow(run_simulation().reshape(grid.shape))\n", - "kernel.kernel.n_threads = 8\n", - "axes[1, 0].imshow(run_simulation().reshape(grid.shape))\n", - "\n", - "axes[1, 1].imshow(run_simulation2());" - ] - } - ], - "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.3" - } - }, - "nbformat": 4, - "nbformat_minor": 5 + { + "cell_type": "markdown", + "id": "076b1ff5-7333-438b-a25f-735434af0543", + "metadata": {}, + "source": ["# Fastscapelib \"user flow kernels\" prototype"] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c79bc75-4c5b-471f-aa79-dd44e1069f1b", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numba as nb\n", + "import numpy as np\n", + "\n", + "import fastscapelib as fs\n", + "\n", + "grid = fs.RasterGrid.from_length([201, 301], [5e4, 7.5e4], fs.NodeStatus.FIXED_VALUE)\n", + "flow_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()])\n", + "\n", + "rng = np.random.Generator(np.random.PCG64(1234))\n", + "\n", + "init_elevation = rng.uniform(0, 5, size=grid.shape)\n", + "drainage_area = np.empty_like(init_elevation)\n", + "\n", + "uplift_rate = np.full_like(init_elevation, 1e-3)\n", + "uplift_rate[[0, -1], :] = 0.0\n", + "uplift_rate[:, [0, -1]] = 0.0\n", + "flow_graph.update_routes(init_elevation)\n", + "\n", + "\n", + "def kernel_func(node): \n", + " dt = node.dt\n", + " r_count = node.receivers.count\n", + " if r_count == 1 and node.receivers.distance[0] == 0.0:\n", + " return\n", + "\n", + " elevation_flooded = np.finfo(np.double).max\n", + "\n", + " for r in range(r_count):\n", + " irec_elevation_next = node.receivers.elevation[r] - node.receivers.erosion[r]\n", + "\n", + " if irec_elevation_next < elevation_flooded:\n", + " elevation_flooded = irec_elevation_next\n", + "\n", + " if node.elevation <= elevation_flooded:\n", + " return\n", + "\n", + " eq_num = node.elevation\n", + " eq_den = 1.0\n", + "\n", + " for r in range(r_count):\n", + " irec_elevation = node.receivers.elevation[r]\n", + " irec_elevation_next = irec_elevation - node.receivers.erosion[r]\n", + "\n", + " if irec_elevation > node.elevation:\n", + " continue\n", + "\n", + " irec_weight = node.receivers.weight[r]\n", + " irec_distance = node.receivers.distance[r]\n", + "\n", + " factor = (\n", + " node.k_coef * node.dt * np.power(node.drainage_area * irec_weight, node.area_exp)\n", + " )\n", + " factor /= irec_distance\n", + " eq_num += factor * irec_elevation_next\n", + " eq_den += factor\n", + "\n", + " elevation_updated = eq_num / eq_den\n", + "\n", + " if elevation_updated < elevation_flooded:\n", + " elevation_updated = elevation_flooded + np.finfo(np.double).tiny\n", + "\n", + " node.erosion = node.elevation - elevation_updated\n", + "\n", + "elevation = init_elevation.ravel().copy()\n", + "erosion = np.zeros(flow_graph.size)\n", + "drainage_area = np.ones(flow_graph.size)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4222ec0-82c0-4617-abef-9dfd44e30ad8", + "metadata": {}, + "outputs": [], + "source": [ + "kernel, data = fs.create_flow_kernel(\n", + " flow_graph,\n", + " kernel_func,\n", + " spec=dict(\n", + " elevation=nb.float64[::1],\n", + " erosion=nb.float64[::1],\n", + " drainage_area=nb.float64[::1],\n", + " k_coef=(nb.float64, 2e-4),\n", + " area_exp=(nb.float64, 0.4),\n", + " slope_exp=(nb.float64, 1.),\n", + " dt=(nb.float64, 2e4),\n", + " ),\n", + " outputs=[\"erosion\"],\n", + " max_receivers=1,\n", + " n_threads=1,\n", + " # print_generated_code=True,\n", + " apply_dir=fs.flow.FlowGraphTraversalDir.BREADTH_UPSTREAM\n", + ")\n", + "\n", + "data.bind(\n", + " elevation=elevation,\n", + " erosion=erosion,\n", + " drainage_area=drainage_area,\n", + ")\n", + "\n", + "flow_graph.apply_kernel(kernel, data)\n", + "\n", + "import time\n", + "N_TIMES = 1\n", + "\n", + "t1 = time.time()\n", + "for _ in range(N_TIMES):\n", + " flow_graph.apply_kernel(kernel, data)\n", + "print(round((time.time() - t1) * 1e9 / N_TIMES / grid.size, 1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17991b9f-69fb-44e3-abed-c46af557a8d1", + "metadata": {}, + "outputs": [], + "source": [ + "data.dt = dt = 2e4\n", + "nsteps = 50\n", + "uplift = dt * uplift_rate.ravel()\n", + "\n", + "def run_simulation():\n", + " elevation = init_elevation.copy().ravel()\n", + " drainage_area = np.empty_like(elevation)\n", + "\n", + " for step in range(nsteps):\n", + " # uplift (no uplift at fixed elevation boundaries)\n", + " uplifted_elevation = elevation + uplift\n", + " \n", + " # flow routing\n", + " flow_graph.update_routes(uplifted_elevation)\n", + "\n", + " # flow accumulation (drainage area)\n", + " flow_graph.accumulate(drainage_area, 1.0)\n", + "\n", + " data.erosion.fill(0.)\n", + " data.bind(elevation=uplifted_elevation, drainage_area=drainage_area, dt=dt)\n", + " flow_graph.apply_kernel(kernel, data)\n", + "\n", + " # update topography\n", + " elevation = uplifted_elevation - data.erosion\n", + "\n", + " return elevation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c80f6143-445e-4d8f-8520-c54b7b86e7bb", + "metadata": {}, + "outputs": [], + "source": [ + "for i in [1, 2, 4, 8]:\n", + " kernel.kernel.n_threads = i\n", + " %timeit -r 10 -n 10 flow_graph.apply_kernel(kernel, data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90de6ce1-4ba7-4225-a60c-6f957b2c2691", + "metadata": {}, + "outputs": [], + "source": [ + "for i in [1, 2, 4, 8]:\n", + " kernel.kernel.n_threads = i\n", + " %timeit -r 1 -n 1 run_simulation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a6598a97-acf5-4a14-a7e6-5310e3688a6d", + "metadata": {}, + "outputs": [], + "source": [ + "# efficiency of the parallel execution depends to the kernel size, the current flow graph and the kernel execution order\n", + "for i in [1, 2, 4, 8]:\n", + " kernel.kernel.n_threads = i\n", + " %timeit -r 10 -n 10 flow_graph.apply_kernel(kernel, data)" + ] + }, + { + "cell_type": "markdown", + "id": "1235bd7e-0dae-465a-ab6f-630b0a6e3cb5", + "metadata": {}, + "source": ["### Fastscapelib Python bindings"] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10c5fae2-dbf7-4971-9c16-73429ff9bca0", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numba as nb\n", + "import numpy as np\n", + "\n", + "import fastscapelib as fs\n", + "\n", + "uplift_rate = np.full_like(init_elevation, 1e-3)\n", + "uplift_rate[[0, -1], :] = 0.\n", + "uplift_rate[:, [0, -1]] = 0.\n", + "\n", + "dt = 2e4\n", + "nsteps = 50\n", + "\n", + "k_coef = 2e-4\n", + "area_exp = 0.4\n", + "slope_exp = 1\n", + "\n", + "flow_graph.update_routes(init_elevation)\n", + "\n", + "spl_eroder = fs.SPLEroder(\n", + " flow_graph,\n", + " k_coef=k_coef,\n", + " area_exp=area_exp,\n", + " slope_exp=slope_exp,\n", + " tolerance=1e-5,\n", + ")\n", + "def run_simulation2():\n", + " elevation = init_elevation.copy()\n", + " drainage_area = np.empty_like(init_elevation)\n", + "\n", + " for step in range(nsteps):\n", + " # uplift (no uplift at fixed elevation boundaries)\n", + " uplifted_elevation = elevation + dt * uplift_rate\n", + " \n", + " # flow routing\n", + " filled_elevation = flow_graph.update_routes(uplifted_elevation)\n", + " \n", + " # flow accumulation (drainage area)\n", + " flow_graph.accumulate(drainage_area, 1.0)\n", + " \n", + " # apply channel erosion (SPL)\n", + " spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, dt)\n", + "\n", + " # update topography\n", + " elevation = uplifted_elevation - spl_erosion\n", + "\n", + " return elevation " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7089e1f-ca26-42a9-89ee-e2711d5216d9", + "metadata": {}, + "outputs": [], + "source": ["%timeit -r 1 -n 1 run_simulation2()"] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8db53d5a-a17d-4079-86d0-bed75bc1d999", + "metadata": {}, + "outputs": [], + "source": ["%timeit -r 10 -n 10 spl_eroder.erode(elevation, drainage_area, dt)"] + }, + { + "cell_type": "markdown", + "id": "a50faaa5-8c71-4d7c-acce-faa733f73c49", + "metadata": {}, + "source": ["### Benchmark Results"] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ab15183-51e6-4072-8e69-a06c24d6757c", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(20, 10))\n", + "\n", + "kernel.kernel.n_threads = 1\n", + "axes[0, 0].imshow(run_simulation().reshape(grid.shape))\n", + "kernel.kernel.n_threads = 4\n", + "axes[0, 1].imshow(run_simulation().reshape(grid.shape))\n", + "kernel.kernel.n_threads = 8\n", + "axes[1, 0].imshow(run_simulation().reshape(grid.shape))\n", + "\n", + "axes[1, 1].imshow(run_simulation2());" + ] + } + ], + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 } diff --git a/include/fastscapelib/flow/flow_router.hpp b/include/fastscapelib/flow/flow_router.hpp index c5783324..73a9d710 100644 --- a/include/fastscapelib/flow/flow_router.hpp +++ b/include/fastscapelib/flow/flow_router.hpp @@ -27,7 +27,7 @@ namespace fastscapelib public: single_flow_router() = default; single_flow_router(int n_threads) - : m_threads_count(n_threads){}; + : m_threads_count(n_threads) {}; inline std::string name() const noexcept override { @@ -70,7 +70,7 @@ namespace fastscapelib public: flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; void apply(graph_impl_type& graph_impl, data_array_type& elevation, @@ -268,7 +268,7 @@ namespace fastscapelib using thread_pool_type = thread_pool; flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; void apply(graph_impl_type& graph_impl, data_array_type& elevation, diff --git a/include/fastscapelib/flow/flow_snapshot.hpp b/include/fastscapelib/flow/flow_snapshot.hpp index fe246b04..3a8884b7 100644 --- a/include/fastscapelib/flow/flow_snapshot.hpp +++ b/include/fastscapelib/flow/flow_snapshot.hpp @@ -75,7 +75,7 @@ namespace fastscapelib using elevation_map = std::map>; flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; void save(const FG& graph_impl, graph_impl_map& graph_impl_snapshots, diff --git a/include/fastscapelib/flow/sink_resolver.hpp b/include/fastscapelib/flow/sink_resolver.hpp index 136af01d..d5538645 100644 --- a/include/fastscapelib/flow/sink_resolver.hpp +++ b/include/fastscapelib/flow/sink_resolver.hpp @@ -65,7 +65,7 @@ namespace fastscapelib using thread_pool_type = thread_pool; flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; void apply(graph_impl_type& graph_impl, data_array_type& elevation, @@ -182,7 +182,7 @@ namespace fastscapelib using thread_pool_type = thread_pool; flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; void apply(graph_impl_type& graph_impl, data_array_type& elevation, diff --git a/include/fastscapelib/grid/base.hpp b/include/fastscapelib/grid/base.hpp index 846f1f77..79950422 100644 --- a/include/fastscapelib/grid/base.hpp +++ b/include/fastscapelib/grid/base.hpp @@ -363,7 +363,7 @@ namespace fastscapelib typename neighbors_cache_type::template storage_type; grid(std::size_t size) - : m_neighbors_indices_cache(neighbors_cache_type(size)){}; + : m_neighbors_indices_cache(neighbors_cache_type(size)) {}; ~grid() = default; const derived_grid_type& derived_grid() const noexcept; diff --git a/python/fastscapelib/flow/numba/flow_kernel.py b/python/fastscapelib/flow/numba/flow_kernel.py index df1cba35..b074f9cd 100644 --- a/python/fastscapelib/flow/numba/flow_kernel.py +++ b/python/fastscapelib/flow/numba/flow_kernel.py @@ -74,11 +74,9 @@ def visit_AugAssign(self, node: ast.AugAssign): class NumbaJittedFunc(Protocol[P, R]): - def __call__(self, *args: P.args, **kwargs: P.kwargs) -> R: - ... + def __call__(self, *args: P.args, **kwargs: P.kwargs) -> R: ... - def get_compile_result(self, sig: Any) -> Any: - ... + def get_compile_result(self, sig: Any) -> Any: ... # simplified type annotation for numba.experimental.jitclass decorated class diff --git a/python/src/flow_graph.cpp b/python/src/flow_graph.cpp index 30e5f3e1..36e58cfb 100644 --- a/python/src/flow_graph.cpp +++ b/python/src/flow_graph.cpp @@ -237,10 +237,10 @@ add_flow_graph_bindings(py::module& m) )doc"); mrouter_op.def_readwrite( "slope_exp", &fs::multi_flow_router::m_slope_exp, "Flow partition slope exponent."); - mrouter_op.def("__repr__", - [](const fs::multi_flow_router& op) { - return "MultiFlowRouter (slope_exp=" + std::to_string(op.m_slope_exp) + ")"; - }); + mrouter_op.def( + "__repr__", + [](const fs::multi_flow_router& op) + { return "MultiFlowRouter (slope_exp=" + std::to_string(op.m_slope_exp) + ")"; }); py::class_; - virtual ~flow_graph_impl_wrapper_base(){}; + virtual ~flow_graph_impl_wrapper_base() {}; virtual bool single_flow() const = 0; @@ -96,7 +96,7 @@ namespace fastscapelib class flow_graph_impl_wrapper : public flow_graph_impl_wrapper_base { public: - virtual ~flow_graph_impl_wrapper(){}; + virtual ~flow_graph_impl_wrapper() {}; flow_graph_impl_wrapper(const std::shared_ptr& graph_impl_ptr) : m_graph_impl_ptr(graph_impl_ptr) @@ -192,7 +192,7 @@ namespace fastscapelib template py_flow_graph_impl(const std::shared_ptr& graph_impl_ptr) : m_wrapper_ptr( - std::make_unique>(graph_impl_ptr)){}; + std::make_unique>(graph_impl_ptr)){}; bool single_flow() const { @@ -388,7 +388,7 @@ namespace fastscapelib using shape_type = data_array_type::shape_type; using data_array_size_type = dynamic_shape_container_t; - virtual ~flow_graph_wrapper_base(){}; + virtual ~flow_graph_wrapper_base() {}; virtual bool single_flow() const = 0; virtual size_type size() const = 0; @@ -442,7 +442,7 @@ namespace fastscapelib = std::make_unique(m_snapshot_graph_ptr->impl_ptr()); } - virtual ~flow_graph_wrapper(){}; + virtual ~flow_graph_wrapper() {}; bool single_flow() const { diff --git a/test/test_flow_operator.cpp b/test/test_flow_operator.cpp index e7ddce30..7929bad7 100644 --- a/test/test_flow_operator.cpp +++ b/test/test_flow_operator.cpp @@ -55,7 +55,7 @@ namespace fastscapelib using base_type = flow_operator_impl_base; flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; }; } @@ -91,7 +91,7 @@ namespace fastscapelib using thread_pool_type = thread_pool; flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; void apply(graph_impl_type& /*graph_impl*/, data_array_type& elevation, diff --git a/test/test_flow_router.cpp b/test/test_flow_router.cpp index c8e541cf..0bbdce72 100644 --- a/test/test_flow_router.cpp +++ b/test/test_flow_router.cpp @@ -45,7 +45,7 @@ namespace fastscapelib using base_type = flow_operator_impl; flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; }; template @@ -56,7 +56,7 @@ namespace fastscapelib using base_type = flow_operator_impl; flow_operator_impl(std::shared_ptr ptr) - : base_type(std::move(ptr)){}; + : base_type(std::move(ptr)) {}; }; } diff --git a/test/test_raster_bishop_grid.cpp b/test/test_raster_bishop_grid.cpp index ead9bea7..144e7ee7 100644 --- a/test/test_raster_bishop_grid.cpp +++ b/test/test_raster_bishop_grid.cpp @@ -105,11 +105,7 @@ namespace fastscapelib /*(r,c)*/ })); } // Bottom-right corner - EXPECT_INDICES(4, - 9, - ({ 38 }), - ({ { 3, 8 } - /*(r,c)*/ })); + EXPECT_INDICES(4, 9, ({ 38 }), ({ { 3, 8 } /*(r,c)*/ })); } } #undef EXPECT_INDICES