Creates UxDataArray.laplacian()

docs/api.rst

@@ -529,6 +529,7 @@ Mathematical Operators
    UxDataArray.divergence
    UxDataArray.gradient
    UxDataArray.integrate
+   UxDataArray.laplacian
    UxDataArray.scalardotgradient
 
 

docs/user-guide/vector_calculus.ipynb

@@ -20,9 +20,10 @@
     "1. Gradient of face-centered scalar fields (zonal and meridional components)\n",
     "2. Curl of vector fields (including curl of a gradient ≈ 0 and a synthetic vortex)\n",
     "3. Divergence of vector fields (including Laplacian as ∇·∇φ, divergence of a vortex ≈ 0, and radial expansion > 0)\n",
-    "4. Scalar dot gradient (advection term **v** · ∇q)\n",
-    "5. Vector calculus identities verification\n",
-    "6. Vector calculus for 2D+ arrays (grid dimension plus other dimension(s) such as \"time\")\n"
+    "4. Laplacian of face-centered scalar fields (∇²q)\n",
+    "5. Scalar dot gradient (advection term **v** · ∇q)\n",
+    "6. Vector calculus identities verification\n",
+    "7. Vector calculus for 2D+ arrays (grid dimension plus other dimension(s) such as \"time\")\n"
    ]
   },
   {
@@ -645,7 +646,9 @@
    "outputs": [],
    "source": [
     "# Compute divergence of the gradient (Laplacian)\n",
-    "# ∇²φ = ∂(∇φ_x)/∂x + ∂(∇φ_y)/∂y\n",
+    "# ∇²g = div(grad(g)) = ∂(∇g_x)/∂x + ∂(∇g_y)/∂y - (∂g/∂y)·tan(φ)/a\n",
+    "# (The final term is the metric term from the divergence operator.)\n",
+    "# (This can also be computed by calling uxds[\"gaussian\"].laplacian() directly.)\n",
     "grad_gauss_u = grad_gauss[\"zonal_gradient\"]\n",
     "grad_gauss_v = grad_gauss[\"meridional_gradient\"]\n",
     "\n",
@@ -758,12 +761,109 @@
     "(radial_vectors + radial_div).opts(shared_axes=False)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "36778b62",
+   "metadata": {},
+   "source": [
+    "## 4. Laplacian (∇²q)\n",
+    "\n",
+    "### Background\n",
+    "\n",
+    "The Laplacian of a scalar field q is the divergence of its gradient, ∇ · ∇q. It measures how much the value of q at a point differs from the average of its surroundings, and it shows up in diffusion, wave, and Poisson equations.\n",
+    "\n",
+    "- **Positive Laplacian**: q at that point is smaller than its surroundings (local minimum)\n",
+    "- **Negative Laplacian**: q at that point is larger than its surroundings (local maximum)\n",
+    "- **Zero Laplacian**: q at that point is balanced against its surroundings (harmonic field)\n",
+    "\n",
+    "On a sphere, with $\\lambda$ longitude, $\\varphi$ latitude, and $a$ the sphere radius (`uxgrid.sphere_radius`),\n",
+    "\n",
+    "$$\n",
+    "\\nabla^2 q = \\frac{\\partial^2 q}{\\partial x^2} + \\frac{\\partial^2 q}{\\partial y^2} - \\frac{\\partial q}{\\partial y} \\frac{\\tan\\varphi}{a}\n",
+    "$$\n",
+    "\n",
+    "The last term comes directly from the metric term in the divergence formula.\n",
+    "\n",
+    "### Usage\n",
+    "\n",
+    "The Laplacian can be computed using the `UxDataArray.laplacian()` method on a face-centered data variable:\n",
+    "\n",
+    "| **Input**           |    **Usage**    | **Output**             |\n",
+    "| ------------------- | :-------------: | ---------------------- |\n",
+    "| Scalar field $q$    | `q.laplacian()` | Scalar field $\\nabla^2 q$ |"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e8e352d8",
+   "metadata": {},
+   "source": [
+    "### Example 1: Laplacian of a Gaussian field"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bdbcb653",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Exactly equivalent to the div_of_gradient from the previous section:\n",
+    "laplacian = uxds[\"gaussian\"].laplacian()\n",
+    "diff = laplacian - div_of_gradient\n",
+    "print(f\"Max difference between laplacian and div_of_gradient: {diff.item().2e}\")\n",
+    "laplacian"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c3beb696",
+   "metadata": {},
+   "source": [
+    "### Example 2: Laplacian of a spherical harmonic\n",
+    "The Laplacian of a spherical harmonic of degree $l$ should theoretically follow the formula:\n",
+    "$$\n",
+    "\\nabla^2 q = -\\frac{l(l+1)}{a^2}\\, q\n",
+    "$$\n",
+    "where $a$ is the sphere radius (`uxgrid.sphere_radius`). Spherical harmonics are eigenfunctions of the Laplacian on the sphere: the result is the original field scaled by $-l(l+1)/a^2$, and it does not depend on the order $m$. This holds for any linear combination of spherical harmonics as well, because the Laplacian is a linear operator.\n",
+    "\n",
+    "$x$, $y$, and $z$ can all be formed from linear combinations of spherical harmonics of degree 1, so taking the Laplacian of them should give results close to 0."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "13cd77d5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "grid = uxds.uxgrid\n",
+    "a = grid.sphere_radius\n",
+    "\n",
+    "# Cartesian x of the face centers (unit sphere), wrapped as a UxDataArray\n",
+    "zero = ux.UxDataArray(xr.DataArray(np.zeros(grid.n_face), dims=\"n_face\"), uxgrid=grid)\n",
+    "x_cart = zero + grid.face_x.values\n",
+    "\n",
+    "# x is a degree-1 combination of spherical harmonics, so laplacian(x) = -2x/a^2\n",
+    "lap_x = x_cart.laplacian()\n",
+    "residual = lap_x - (-2.0 * x_cart / a**2)\n",
+    "\n",
+    "finite = np.isfinite(residual.values)\n",
+    "res_vals = residual.values[finite]\n",
+    "print(f\"Residual range: [{res_vals.min():.2e}, {res_vals.max():.2e}]\")\n",
+    "print(f\"Max |residual| / max |analytic|: {np.abs(res_vals).max() / (2 * np.abs(x_cart.values[finite]).max() / a**2):.2%}\")\n",
+    "\n",
+    "residual.plot(cmap=\"RdBu_r\", aspect=1).opts(\n",
+    "    title=\"Laplacian of x minus (-2x/a²) (Should ≈ 0)\", colorbar=True\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "26b6a175",
    "metadata": {},
    "source": [
-    "## 4. Scalar Dot Gradient (**v** · ∇q)\n",
+    "## 5. Scalar Dot Gradient (**v** · ∇q)\n",
     "\n",
     "### Background\n",
     "\n",
@@ -833,7 +933,7 @@
    "id": "35",
    "metadata": {},
    "source": [
-    "## 5. Vector Calculus Identities\n",
+    "## 6. Vector Calculus Identities\n",
     "\n",
     "Let's verify some fundamental vector calculus identities using our computed fields:\n",
     "\n",
@@ -906,7 +1006,7 @@
    "id": "9d247e97",
    "metadata": {},
    "source": [
-    "## 6. Vector calculus for 2D+ arrays\n",
+    "## 7. Vector calculus for 2D+ arrays\n",
     "\n",
     "Real data often has multiple dimensions such as \"time\" or \"levels\", not just the grid dimension.\n",
     "\n",