Comparison with PyPDE: Integrate Benchmarking and Compatibility Analysis

.github/workflows/pdoc.yml

@@ -25,7 +25,7 @@ jobs:
           persist-credentials: false
       - uses: actions/setup-python@v5.2.0
         with:
-          python-version: 3.9
+          python-version: "3.13"
       - env:
           JUPYTER_PLATFORM_DIRS: 1
         run: |

.github/workflows/tests.yml

@@ -38,7 +38,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v5.2.0
       with:
-        python-version: 3.11
+        python-version: "3.13"
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip

examples/PyMPDATA_examples/comparison_against_pypde_2025/cartesian_grid-1.ipynb

@@ -0,0 +1,256 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "[![preview notebook](https://img.shields.io/static/v1?label=render%20on&logo=github&color=87ce3e&message=GitHub)](https://github.com/open-atmos/PyMPDATA/blob/main/examples/PyMPDATA_examples/comparison_against_pypde_2025/cartesian_grid-1.ipynb)\n",
+    "[![launch on Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/drive/1an38vReLVS7q1Wd-qE_CiDtjmPJK-lU1?usp=sharing)\n",
+    " "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "license: GPL v3   \n",
+    "authors: Piotr Karaś, Norbert Klockiewicz, Jakub Kot, Kacper Majchrzak   \n",
+    "Diffusion on a Cartesian grid   \n",
+    "This example shows how to solve the diffusion equation on a Cartesian grid using both py-pde and PyMPDATA."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "if 'google.colab' in sys.modules:\n",
+    "    !pip --quiet install open-atmos-jupyter-utils\n",
+    "    from open_atmos_jupyter_utils import pip_install_on_colab\n",
+    "    pip_install_on_colab('PyMPDATA-examples')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from examples.PyMPDATA_examples.comparison_against_pypde_2025.diffusion_2d import InitialConditions\n",
+    "\n",
+    "initial_conditions = InitialConditions(\n",
+    "    diffusion_coefficient=0.1,\n",
+    "    time_step=0.001,\n",
+    "    time_end=1,\n",
+    "    grid_shape=(30, 18),\n",
+    "    grid_range_x=(-1.0, 1.0),\n",
+    "    grid_range_y=(0.0, 2.0),\n",
+    "    pulse_position=(0.0, 1.0),\n",
+    "    pulse_shape=(2, 2),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2025-06-20T18:48:41.195420Z",
+     "start_time": "2025-06-20T18:48:41.055293Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from examples.PyMPDATA_examples.comparison_against_pypde_2025.diffusion_2d import create_pde_state\n",
+    "\n",
+    "from open_atmos_jupyter_utils import show_plot\n",
+    "from pde import DiffusionPDE\n",
+    "\n",
+    "state = create_pde_state(initial_conditions=initial_conditions)\n",
+    "state.plot(title=\"Initial condition\", cmap=\"magma\")\n",
+    "show_plot(filename=\"Initial_condition\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "eq = DiffusionPDE(0.1)  # define the pde\n",
+    "result_pdepde = eq.solve(state, t_range=1, dt=0.001)\n",
+    "result_pdepde.plot(cmap=\"magma\")\n",
+    "result_pdepde = result_pdepde.data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2025-06-20T18:48:45.524859Z",
+     "start_time": "2025-06-20T18:48:44.260427Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from examples.PyMPDATA_examples.comparison_against_pypde_2025.diffusion_2d import py_pde_solution\n",
+    "\n",
+    "\n",
+    "result_pdepde2 = py_pde_solution(initial_conditions=initial_conditions)\n",
+    "assert (result_pdepde == result_pdepde2).all()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2025-06-20T18:50:59.449433Z",
+     "start_time": "2025-06-20T18:50:58.089985Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt\n",
+    "from examples.PyMPDATA_examples.comparison_against_pypde_2025.diffusion_2d import (\n",
+    "    mpdata_solution,\n",
+    ")\n",
+    "\n",
+    "result_mpdata = mpdata_solution(initial_conditions=initial_conditions)\n",
+    "\n",
+    "# Plot final result\n",
+    "plt.figure(figsize=(6, 4))\n",
+    "plt.imshow(\n",
+    "    result_mpdata,\n",
+    "    cmap=\"magma\",\n",
+    "    origin=\"lower\",\n",
+    "    extent=[\n",
+    "        initial_conditions.min_x,\n",
+    "        initial_conditions.max_x,\n",
+    "        initial_conditions.min_y,\n",
+    "        initial_conditions.max_y,\n",
+    "    ],\n",
+    ")\n",
+    "plt.title(\"PyMPDATA Diffusion Result\")\n",
+    "plt.xlabel(\"x\")\n",
+    "plt.ylabel(\"y\")\n",
+    "plt.colorbar()\n",
+    "plt.tight_layout()\n",
+    "show_plot(filename=\"PyMPDATA_diffusion_result\")\n",
+    "\n",
+    "# Total mass check\n",
+    "mass = np.sum(result_mpdata) * initial_conditions.dx * initial_conditions.dy\n",
+    "print(f\"Total mass after diffusion: {mass:.6f} (should be close to 1)\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2025-06-20T18:51:09.076355Z",
+     "start_time": "2025-06-20T18:51:09.072679Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "diff = result_mpdata - result_pdepde\n",
+    "rmse = np.sqrt(np.mean(diff**2))\n",
+    "l1 = np.sum(np.abs(diff)) * initial_conditions.dx * initial_conditions.dy\n",
+    "mass_py_pde = result_pdepde.sum() * initial_conditions.dx * initial_conditions.dy\n",
+    "mass_mpdata = result_mpdata.sum() * initial_conditions.dx * initial_conditions.dy\n",
+    "print(\"Total mass:\")\n",
+    "print(f\"py-pde:{mass_py_pde:>12.6f}\")\n",
+    "print(f\"PyMPDATA:{mass_mpdata:>10.6f}\")\n",
+    "\n",
+    "print(f\"RMSE: {rmse:.6f}, L1 Error: {l1:.6f}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2025-06-20T18:53:49.719836Z",
+     "start_time": "2025-06-20T18:53:49.672914Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def analytical_solution(x, y, t, D=0.1, center=(0.0, 1.0)):\n",
+    "    x0, y0 = center\n",
+    "    r2 = (x - x0) ** 2 + (y - y0) ** 2\n",
+    "    return (1 / (4 * np.pi * D * t)) * np.exp(-r2 / (4 * D * t))\n",
+    "\n",
+    "\n",
+    "x = np.linspace(initial_conditions.min_x + initial_conditions.dx / 2, initial_conditions.max_x - initial_conditions.dx / 2, initial_conditions.nx)\n",
+    "y = np.linspace(initial_conditions.min_y + initial_conditions.dy / 2, initial_conditions.max_y - initial_conditions.dy / 2, initial_conditions.ny)\n",
+    "X, Y = np.meshgrid(x, y, indexing=\"ij\")\n",
+    "u_exact = analytical_solution(X, Y, initial_conditions.time_end, D=initial_conditions.diffusion_coefficient)\n",
+    "\n",
+    "plt.figure(figsize=(6, 4))\n",
+    "plt.imshow(u_exact.T, origin=\"lower\", extent=[\n",
+    "    initial_conditions.min_x,\n",
+    "    initial_conditions.max_x,\n",
+    "    initial_conditions.min_y,\n",
+    "    initial_conditions.max_y,\n",
+    "], cmap=\"magma\")\n",
+    "plt.title(\"Analytical Diffusion Solution\")\n",
+    "plt.colorbar()\n",
+    "show_plot(filename=\"Analytical_diffusion_solution\")\n",
+    "\n",
+    "diff_pde = result_pdepde - u_exact\n",
+    "rmse_pde = np.sqrt(np.mean((diff_pde) ** 2))\n",
+    "l1_error_pde = np.sum(np.abs(diff_pde)) * initial_conditions.dx * initial_conditions.dy\n",
+    "\n",
+    "COL_WIDTH_METHOD = 10\n",
+    "COL_WIDTH_METRIC = 18\n",
+    "\n",
+    "print(\"=\" * 50)\n",
+    "print(\"Comparison with Analytical Solution\")\n",
+    "print(\"=\" * 50)\n",
+    "print(\n",
+    "    f\"{'Method':<{COL_WIDTH_METHOD}} {'RMSE':>{COL_WIDTH_METRIC}} {'L1 Error':>{COL_WIDTH_METRIC}}\"\n",
+    ")\n",
+    "print(\"-\" * 50)\n",
+    "print(\n",
+    "    f\"{'PDE':<{COL_WIDTH_METHOD}} {rmse_pde:>{COL_WIDTH_METRIC}.6f} {l1_error_pde:>{COL_WIDTH_METRIC}.6f}\"\n",
+    ")\n",
+    "\n",
+    "diff_mpdata = result_mpdata - u_exact\n",
+    "rmse_mpdata = np.sqrt(np.mean((diff_mpdata) ** 2))\n",
+    "l1_error_mpdata = np.sum(np.abs(diff_mpdata)) * initial_conditions.dx * initial_conditions.dy\n",
+    "\n",
+    "print(\n",
+    "    f\"{'MPDATA':<{COL_WIDTH_METHOD}} {rmse_mpdata:>{COL_WIDTH_METRIC}.6f} {l1_error_mpdata:>{COL_WIDTH_METRIC}.6f}\"\n",
+    ")\n",
+    "print(\"=\" * 50)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": ".venv",
+   "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.11.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}