From fb6b7b99689b8ea18ca43f154f99f2abe04ce74f Mon Sep 17 00:00:00 2001 From: Mark Piper Date: Mon, 31 Aug 2026 15:41:58 -0600 Subject: [PATCH] Add Landlab diffusion model and example notebooks --- extra/landlab_diffusion_model.ipynb | 154 ++++++++++++++++++++++++ extra/landlab_diffusion_model.py | 51 ++++++++ extra/run_landlab_diffusion_model.ipynb | 53 ++++++++ 3 files changed, 258 insertions(+) create mode 100644 extra/landlab_diffusion_model.ipynb create mode 100644 extra/landlab_diffusion_model.py create mode 100644 extra/run_landlab_diffusion_model.ipynb diff --git a/extra/landlab_diffusion_model.ipynb b/extra/landlab_diffusion_model.ipynb new file mode 100644 index 0000000..a9f729d --- /dev/null +++ b/extra/landlab_diffusion_model.ipynb @@ -0,0 +1,154 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "12022014-5f95-4708-8d40-c16f13353b91", + "metadata": {}, + "source": [ + "# Write a 2D diffusion model with Landlab" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "83ecd54c-376b-4a80-8c59-7a6ff1029afa", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "import tomllib\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6aef160f-bf7f-494b-ba49-ba61f2f06a9f", + "metadata": {}, + "outputs": [], + "source": [ + "from landlab import RasterModelGrid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "442c2d18-0a1d-43ee-9951-dfb49a49f47c", + "metadata": {}, + "outputs": [], + "source": [ + "def calculate_stable_time_step(dx, diffusivity):\n", + " return 0.25 * dx**2 / diffusivity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4122db3c-2f16-4783-aa70-e0ff439cef68", + "metadata": {}, + "outputs": [], + "source": [ + "def new_profile(grid, step_at=1.0):\n", + " z = grid.zeros(at=\"node\")\n", + " z[grid.x_of_node >= step_at] = 1.0\n", + " return z" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7513667-d45e-4dee-bf1c-7d3a315db973", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_profile(grid, concentration, color=\"r\"):\n", + " grid.imshow(concentration)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a3b0cdbf-9c56-44db-b1b1-57b1e950eb76", + "metadata": {}, + "outputs": [], + "source": [ + "def calculate_elevation_change(grid, z, diffusivity):\n", + " dzdl = grid.calc_grad_at_link(z)\n", + " qs_at_link = -diffusivity * dzdl\n", + " dzdt = -grid.calc_flux_div_at_node(qs_at_link)\n", + " return dzdt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa1d8b0d-25e3-46ed-8c82-24b24c3ba8c1", + "metadata": {}, + "outputs": [], + "source": [ + "def diffuse_until(grid, z_initial, stop_time, diffusivity=1.0):\n", + " stable_dt = 0.9 * calculate_stable_time_step(np.min(grid.length_of_link), diffusivity)\n", + " z = z_initial.copy()\n", + " \n", + " time = 0\n", + " while time < stop_time:\n", + " dt = min(stable_dt, stop_time - time)\n", + " z += dt* calculate_elevation_change(grid, z, diffusivity)\n", + " time += dt\n", + "\n", + " return z" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "623c29b7-6d26-41ce-8041-58fbaa0784dd", + "metadata": {}, + "outputs": [], + "source": [ + "def run_diffusion_model():\n", + " shape = (100, 200)\n", + " stop_time = 5.0\n", + " diffusivity = 10.0\n", + " grid = RasterModelGrid(shape, xy_spacing=(1.0, 1.0))\n", + " z_initial = new_profile(grid, step_at=100)\n", + " z = diffuse_until(grid, z_initial, stop_time, diffusivity=diffusivity)\n", + " plot_profile(grid, z)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6cfcd5c-9660-429f-9705-eadc2076b612", + "metadata": {}, + "outputs": [], + "source": [ + "run_diffusion_model()" + ] + } + ], + "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.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/extra/landlab_diffusion_model.py b/extra/landlab_diffusion_model.py new file mode 100644 index 0000000..3d0eba9 --- /dev/null +++ b/extra/landlab_diffusion_model.py @@ -0,0 +1,51 @@ +import os +import sys +import tomllib + +import numpy as np +import matplotlib.pyplot as plt +from landlab import RasterModelGrid + + +def calculate_stable_time_step(dx, diffusivity): + return 0.25 * dx**2 / diffusivity + + +def new_profile(grid, step_at=1.0): + z = grid.zeros(at="node") + z[grid.x_of_node >= step_at] = 1.0 + return z + + +def plot_profile(grid, concentration, color="r"): + grid.imshow(concentration) + + +def calculate_elevation_change(grid, z, diffusivity): + dzdl = grid.calc_grad_at_link(z) + qs_at_link = -diffusivity * dzdl + dzdt = -grid.calc_flux_div_at_node(qs_at_link) + return dzdt + + +def diffuse_until(grid, z_initial, stop_time, diffusivity=1.0): + stable_dt = 0.9 * calculate_stable_time_step(np.min(grid.length_of_link), diffusivity) + z = z_initial.copy() + + time = 0 + while time < stop_time: + dt = min(stable_dt, stop_time - time) + z += dt* calculate_elevation_change(grid, z, diffusivity) + time += dt + + return z + + +def run_diffusion_model(): + shape = (100, 200) + stop_time = 5.0 + diffusivity = 10.0 + grid = RasterModelGrid(shape, xy_spacing=(1.0, 1.0)) + z_initial = new_profile(grid, step_at=100) + z = diffuse_until(grid, z_initial, stop_time, diffusivity=diffusivity) + plot_profile(grid, z) diff --git a/extra/run_landlab_diffusion_model.ipynb b/extra/run_landlab_diffusion_model.ipynb new file mode 100644 index 0000000..37b8465 --- /dev/null +++ b/extra/run_landlab_diffusion_model.ipynb @@ -0,0 +1,53 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1e9709b6-0a74-43cd-8d7a-3116200b6105", + "metadata": {}, + "source": [ + "# Run the Landlab diffusion model from a script" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48f8c130-b687-4f6c-88ea-3a26888cc06b", + "metadata": {}, + "outputs": [], + "source": [ + "from landlab_diffusion_model import run_diffusion_model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9eb8d8f-f1c6-459e-9e4f-55e14acd2b64", + "metadata": {}, + "outputs": [], + "source": [ + "run_diffusion_model()" + ] + } + ], + "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.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}