Add Landlab diffusion model and example notebooks

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
+}

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)

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
+}