From 609f80e0231ee8d77f007b0d8d2063ded5d5fe21 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 3 Sep 2026 10:55:21 -0500 Subject: [PATCH 1/6] Keep divergence lazy on dask-backed input The spherical metric term added in #1663 calls `other.values`, which materializes the field and forces the whole divergence to compute eagerly. Its companion in `curl` uses `.data` and stays lazy, so switch divergence to match. --- uxarray/core/dataarray.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index 1343d07f6..0220175f9 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -1796,7 +1796,7 @@ def divergence( # Spherical metric term, the companion of the one in curl(). Omitting # it is only valid on a plane. tan_lat = np.tan(np.deg2rad(self.uxgrid.face_lat.values)) - metric = other.values * tan_lat + metric = other.data * tan_lat if scale_by_radius and "sphere_radius" in self.uxgrid._ds.attrs: metric = metric / self.uxgrid._ds.attrs["sphere_radius"] divergence = divergence - metric From 248bb8c95a77faf254cc029aada1f57e1fe687c7 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 3 Sep 2026 11:20:25 -0500 Subject: [PATCH 2/6] Teach the spherical operators in the vector calculus user guide The curl and divergence background cells give the planar formulas, and eight code cells hand-roll the operators out of gradient() components instead of calling curl()/divergence(). Since #1663 the library adds the spherical metric terms, so every worked example demonstrated something the library no longer computes. Add the metric terms to the math, switch the examples to the public API, and rewrite the constant-field discussion: curl and divergence of a constant field are u*tan(lat)/a and -v*tan(lat)/a, not zero. Both are now checked against their closed forms. Residual magnitudes in the prose are corrected to match what the notebook actually prints. --- docs/user-guide/vector_calculus.ipynb | 129 ++++++++++---------------- 1 file changed, 48 insertions(+), 81 deletions(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index 27ef2ed39..3cce25140 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/docs/user-guide/vector_calculus.ipynb @@ -126,7 +126,7 @@ "id": "66caacb5", "metadata": {}, "source": [ - "> **Units.** By default, `gradient()` and `curl()` divide by `uxgrid.sphere_radius` so derivatives carry physical units (e.g. `[data units]/m`, `1/s` for velocity curl). Pass `scale_by_radius=False` to keep results on the unit sphere (per radian). If the grid has no `sphere_radius` attribute, the call falls back to unit-sphere output and emits a `UserWarning`." + "> **Units.** By default, `gradient()`, `curl()`, and `divergence()` divide by `uxgrid.sphere_radius` so derivatives carry physical units (e.g. `[data units]/m`, `1/s` for velocity curl). Pass `scale_by_radius=False` to keep results on the unit sphere (per radian). If the grid has no `sphere_radius` attribute, the call falls back to unit-sphere output and emits a `UserWarning`." ] }, { @@ -228,14 +228,16 @@ "\n", "### Background\n", "\n", - "The curl of a vector field **F** = (u, v) measures the local rotation or circulation. In 2D, curl produces a scalar field representing the magnitude of rotation:\n", + "The curl of a vector field **F** = (u, v) measures the local rotation or circulation. In 2D it produces a scalar field. On a sphere of radius $a$ it carries a metric term that the planar formula leaves out:\n", "\n", "$$\n", - "\\text{curl}(\\mathbf{F}) = \\nabla \\times \\mathbf{F} = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y}\n", + "\\text{curl}(\\mathbf{F}) = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y} + \\frac{u \\tan\\varphi}{a}\n", "$$\n", "\n", + "That last term comes from the convergence of the meridians and vanishes only on the equator. Dropping it costs a factor of two on solid-body rotation, so `curl()` includes it. Building the same quantity by hand out of `gradient()` components gives the planar answer instead.\n", + "\n", "- **Positive curl**: Counter-clockwise rotation\n", - "- **Negative curl**: Clockwise rotation \n", + "- **Negative curl**: Clockwise rotation\n", "- **Zero curl**: No local rotation (irrotational flow)\n", "\n", "### Usage\n", @@ -244,7 +246,7 @@ "\n", "| **Input** | **Usage** | **Output** |\n", "| ---------------------------- | :---------------------: | -------------------- |\n", - "| Vector field (u, v) | `u.curl(v)` | Scalar curl field |" + "| Vector field (u, v) | `u.curl(v)` | Scalar curl field |\n" ] }, { @@ -254,7 +256,9 @@ "source": [ "### Constant Fields (Mathematical Validation)\n", "\n", - "The curl of a constant vector field should be zero everywhere (within numerical precision)." + "On a plane the curl of a constant vector field is zero everywhere. On a sphere it is not: with $u$ constant, both derivative terms drop out and the metric term survives, leaving $\\zeta = u\\tan\\varphi/a$ exactly.\n", + "\n", + "This grid subset spans about ±2° of latitude, so $\\tan\\varphi/a$ stays below $5.5\\times10^{-9}\\ \\mathrm{m^{-1}}$ and the result looks like round-off — but it is a real signal, and we can check it against the closed form.\n" ] }, { @@ -268,12 +272,7 @@ "u_constant = uxds[\"face_lat\"] * 0 + 1.0\n", "v_constant = uxds[\"face_lat\"] * 0 + 2.0\n", "\n", - "# Compute partials via gradient\n", - "grad_u = u_constant.gradient()\n", - "grad_v = v_constant.gradient()\n", - "du_dy = grad_u[\"meridional_gradient\"]\n", - "dv_dx = grad_v[\"zonal_gradient\"]\n", - "curl_constant = dv_dx - du_dy\n", + "curl_constant = u_constant.curl(v_constant)\n", "\n", "finite = np.isfinite(curl_constant.values)\n", "vals = curl_constant.values[finite]\n", @@ -284,7 +283,12 @@ " print(f\"Finite curl range: [{vals.min():.2e}, {vals.max():.2e}]\")\n", " print(\n", " f\"Max |curl|: {np.abs(vals).max():.2e}, Mean |curl|: {np.abs(vals).mean():.2e}\"\n", - " )" + " )\n", + "\n", + "# The closed form for a constant u on the sphere is u*tan(lat)/a\n", + "lat_rad = np.deg2rad(uxds.uxgrid.face_lat.values)\n", + "curl_exact = 1.0 * np.tan(lat_rad) / uxds.uxgrid.sphere_radius\n", + "print(f\"Max deviation from u*tan(lat)/a: {np.abs(vals - curl_exact[finite]).max():.2e}\")" ] }, { @@ -306,12 +310,7 @@ "u_gauss = uxds[\"gaussian\"]\n", "v_gauss = uxds[\"inverse_gaussian\"]\n", "\n", - "# Compute partials via gradient\n", - "grad_u = u_gauss.gradient()\n", - "grad_v = v_gauss.gradient()\n", - "du_dy = grad_u[\"meridional_gradient\"]\n", - "dv_dx = grad_v[\"zonal_gradient\"]\n", - "curl_gauss = dv_dx - du_dy\n", + "curl_gauss = u_gauss.curl(v_gauss)\n", "\n", "finite = np.isfinite(curl_gauss.values)\n", "vals = curl_gauss.values[finite]\n", @@ -319,8 +318,8 @@ " f\"Total faces: {curl_gauss.size}, interior: {vals.size}, boundary NaNs: {np.isnan(curl_gauss.values).sum()}\"\n", ")\n", "if vals.size:\n", - " print(f\"Finite curl range: [{vals.min():.6f}, {vals.max():.6f}]\")\n", - " print(f\"Mean curl (finite): {vals.mean():.6f}\")" + " print(f\"Finite curl range: [{vals.min():.6e}, {vals.max():.6e}]\")\n", + " print(f\"Mean curl (finite): {vals.mean():.6e}\")" ] }, { @@ -332,7 +331,7 @@ "\n", "In the continuous setting, curl(∇φ) = 0 exactly for any scalar field φ. On an unstructured mesh, however, gradient and curl are independent finite-volume stencils that do not form a discrete de Rham complex, so the identity holds only approximately. The residual is a **discretization error** — not a bug — and shrinks with grid refinement.\n", "\n", - "The *magnitude* of the residual depends on the units. By default `gradient()`/`curl()` scale by `uxgrid.sphere_radius` (Earth ≈ 6.37×10⁶ m), so each derivative picks up a factor of 1/radius. `curl(∇φ)` applies the gradient stencil twice and therefore carries a factor of 1/radius² (≈ 4×10⁻¹⁴). For the Gaussian field here (φ ~ O(1)) the scaled residual is ~O(10⁻¹³); on the unit sphere (`scale_by_radius=False`) the same residual is ~O(1–10). Either way it shrinks with refinement.\n" + "The *magnitude* of the residual depends on the units. By default `gradient()`/`curl()` scale by `uxgrid.sphere_radius` (Earth ≈ 6.37×10⁶ m), so each derivative picks up a factor of 1/radius. `curl(∇φ)` applies the gradient stencil twice and therefore carries a factor of 1/radius² (≈ 2×10⁻¹⁴). For the Gaussian field here (φ ~ O(1)) the scaled residual is ~O(10⁻¹⁵); on the unit sphere (`scale_by_radius=False`) the same residual is ~O(0.1). Either way it shrinks with refinement.\n" ] }, { @@ -342,25 +341,18 @@ "metadata": {}, "outputs": [], "source": [ - "# Extract gradient components\n", + "# Extract gradient components and treat them as a vector field\n", "u_component = grad_gauss.zonal_gradient\n", "v_component = grad_gauss.meridional_gradient\n", "\n", - "# Compute partial derivatives via gradient()\n", - "grad_u = u_component.gradient()\n", - "grad_v = v_component.gradient()\n", - "du_dy = grad_u[\"meridional_gradient\"]\n", - "dv_dx = grad_v[\"zonal_gradient\"]\n", - "\n", - "# Curl = ∂v/∂x - ∂u/∂y\n", - "curl_of_gradient = dv_dx - du_dy\n", + "curl_of_gradient = u_component.curl(v_component)\n", "\n", "print(\n", " f\"Curl of gradient range: [{curl_of_gradient.min().values:.2e}, {curl_of_gradient.max().values:.2e}]\"\n", ")\n", "print(f\"Mean absolute curl: {abs(curl_of_gradient).mean().values:.2e}\")\n", "\n", - "# Note: values are ~O(1e-13) (per meter^2) with the default radius scaling,\n", + "# Note: values are ~O(1e-15) (per meter^2) with the default radius scaling,\n", "# so we let the color limits autoscale rather than hardcoding them.\n", "curl_plot = curl_of_gradient.plot(cmap=\"RdBu_r\", aspect=1).opts(\n", " title=\"Curl of Gradient Field (Should ≈ 0)\", colorbar=True\n", @@ -408,12 +400,8 @@ " v_vortex_data, dims=[\"n_face\"], uxgrid=uxds.uxgrid, name=\"v_vortex\"\n", ")\n", "\n", - "# Compute curl via gradients (default: scaled by sphere_radius -> per meter)\n", - "grad_u = u_vortex.gradient()\n", - "grad_v = v_vortex.gradient()\n", - "du_dy = grad_u[\"meridional_gradient\"]\n", - "dv_dx = grad_v[\"zonal_gradient\"]\n", - "curl_vortex = dv_dx - du_dy\n", + "# Compute curl (default: scaled by sphere_radius -> per meter)\n", + "curl_vortex = u_vortex.curl(v_vortex)\n", "\n", "print(\n", " f\"Vortex curl range: [{curl_vortex.min().values:.2e}, {curl_vortex.max().values:.2e}]\"\n", @@ -546,10 +534,10 @@ "\n", "### Background\n", "\n", - "The divergence of a vector field **F** = (u, v) measures the local expansion or contraction of the field:\n", + "The divergence of a vector field **F** = (u, v) measures the local expansion or contraction of the field. It carries the companion of the metric term in `curl`, with the opposite sign and the meridional component in place of the zonal one:\n", "\n", "$$\n", - "\\text{div}(\\mathbf{F}) = \\nabla \\cdot \\mathbf{F} = \\frac{\\partial u}{\\partial x} + \\frac{\\partial v}{\\partial y}\n", + "\\text{div}(\\mathbf{F}) = \\frac{\\partial u}{\\partial x} + \\frac{\\partial v}{\\partial y} - \\frac{v \\tan\\varphi}{a}\n", "$$\n", "\n", "- **Positive divergence**: Expansion (source)\n", @@ -562,7 +550,7 @@ "\n", "| **Input** | **Usage** | **Output** |\n", "| ---------------------------- | :---------------------: | ----------------------- |\n", - "| Vector field (u, v) | `u.divergence(v)` | Scalar divergence field |" + "| Vector field (u, v) | `u.divergence(v)` | Scalar divergence field |\n" ] }, { @@ -572,7 +560,7 @@ "source": [ "### Constant Fields (Mathematical Validation)\n", "\n", - "The divergence of a constant vector field should be zero everywhere (within numerical precision)." + "As with curl, a constant field is divergence-free on a plane but not on a sphere. The derivative terms vanish and the metric term leaves $-v\\tan\\varphi/a$, which is twice the curl residual above and of the opposite sign, since $v = 2u$ here.\n" ] }, { @@ -586,12 +574,7 @@ "u_constant = uxds[\"face_lat\"] * 0 + 1.0\n", "v_constant = uxds[\"face_lat\"] * 0 + 2.0\n", "\n", - "# Compute partials via gradient\n", - "grad_u = u_constant.gradient()\n", - "grad_v = v_constant.gradient()\n", - "du_dx = grad_u[\"zonal_gradient\"]\n", - "dv_dy = grad_v[\"meridional_gradient\"]\n", - "div_constant = du_dx + dv_dy\n", + "div_constant = u_constant.divergence(v_constant)\n", "\n", "finite = np.isfinite(div_constant.values)\n", "vals = div_constant.values[finite]\n", @@ -600,7 +583,11 @@ ")\n", "if vals.size:\n", " print(f\"Finite divergence range: [{vals.min():.2e}, {vals.max():.2e}]\")\n", - " print(f\"Max |div|: {np.abs(vals).max():.2e}, Mean |div|: {np.abs(vals).mean():.2e}\")" + " print(f\"Max |div|: {np.abs(vals).max():.2e}, Mean |div|: {np.abs(vals).mean():.2e}\")\n", + "\n", + "# The closed form for a constant v on the sphere is -v*tan(lat)/a\n", + "div_exact = -2.0 * np.tan(lat_rad) / uxds.uxgrid.sphere_radius\n", + "print(f\"Max deviation from -v*tan(lat)/a: {np.abs(vals - div_exact[finite]).max():.2e}\")" ] }, { @@ -622,12 +609,7 @@ "u_gauss = uxds[\"gaussian\"]\n", "v_gauss = uxds[\"inverse_gaussian\"]\n", "\n", - "# Compute partials via gradient\n", - "grad_u = u_gauss.gradient()\n", - "grad_v = v_gauss.gradient()\n", - "du_dx = grad_u[\"zonal_gradient\"]\n", - "dv_dy = grad_v[\"meridional_gradient\"]\n", - "div_gauss = du_dx + dv_dy\n", + "div_gauss = u_gauss.divergence(v_gauss)\n", "\n", "finite = np.isfinite(div_gauss.values)\n", "vals = div_gauss.values[finite]\n", @@ -635,8 +617,8 @@ " f\"Total faces: {div_gauss.size}, interior: {vals.size}, boundary NaNs: {np.isnan(div_gauss.values).sum()}\"\n", ")\n", "if vals.size:\n", - " print(f\"Finite divergence range: [{vals.min():.6f}, {vals.max():.6f}]\")\n", - " print(f\"Mean divergence (finite): {vals.mean():.6f}\")" + " print(f\"Finite divergence range: [{vals.min():.6e}, {vals.max():.6e}]\")\n", + " print(f\"Mean divergence (finite): {vals.mean():.6e}\")" ] }, { @@ -654,15 +636,12 @@ "metadata": {}, "outputs": [], "source": [ - "# Compute divergence of the gradient (Laplacian) via partials\n", + "# Compute divergence of the gradient (Laplacian)\n", "# ∇²φ = ∂(∇φ_x)/∂x + ∂(∇φ_y)/∂y\n", "grad_gauss_u = grad_gauss[\"zonal_gradient\"]\n", "grad_gauss_v = grad_gauss[\"meridional_gradient\"]\n", "\n", - "gxu = grad_gauss_u.gradient()[\"zonal_gradient\"] # ∂u/∂x\n", - "gyv = grad_gauss_v.gradient()[\"meridional_gradient\"] # ∂v/∂y\n", - "\n", - "div_of_gradient = gxu + gyv\n", + "div_of_gradient = grad_gauss_u.divergence(grad_gauss_v)\n", "\n", "print(\n", " f\"Divergence of gradient range: [{div_of_gradient.min().values:.2e}, {div_of_gradient.max().values:.2e}]\"\n", @@ -682,7 +661,7 @@ "source": [ "### Example 2: Divergence of Vortex Field\n", "\n", - "Pure rotation should have zero divergence (incompressible):" + "Pure rotation is incompressible, so the divergence should be ≈ 0. The metric term is not identically zero here — $v$ varies across the subset — so the residual sits an order of magnitude above the pure discretization error of the planar formula, and is still negligible against the ~10⁻⁵ signal in the radial case below.\n" ] }, { @@ -692,12 +671,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Compute divergence of the vortex via gradients: div = ∂u/∂x + ∂v/∂y\n", - "grad_u = u_vortex.gradient()\n", - "grad_v = v_vortex.gradient()\n", - "du_dx = grad_u[\"zonal_gradient\"]\n", - "dv_dy = grad_v[\"meridional_gradient\"]\n", - "div_vortex = du_dx + dv_dy\n", + "div_vortex = u_vortex.divergence(v_vortex)\n", "\n", "print(\n", " f\"Vortex divergence range: [{div_vortex.min().values:.2e}, {div_vortex.max().values:.2e}]\"\n", @@ -705,7 +679,7 @@ "print(f\"Mean absolute divergence: {abs(div_vortex).mean().values:.2e}\")\n", "print(\"Pure rotation should have zero divergence\")\n", "\n", - "# Residual is ~O(1e-13) with default radius scaling; let color limits autoscale.\n", + "# Residual is ~O(1e-9) with default radius scaling; let color limits autoscale.\n", "div_vortex_plot = div_vortex.plot(cmap=\"RdBu_r\", aspect=1).opts(\n", " title=\"Divergence of Vortex (Should ≈ 0)\", colorbar=True\n", ")\n", @@ -740,16 +714,9 @@ " v_radial_data, dims=[\"n_face\"], uxgrid=uxds.uxgrid, name=\"v_radial\"\n", ")\n", "\n", - "# Compute curl and divergence via gradients (default: scaled by sphere_radius)\n", - "grad_u = u_radial.gradient()\n", - "grad_v = v_radial.gradient()\n", - "du_dy = grad_u[\"meridional_gradient\"]\n", - "dv_dx = grad_v[\"zonal_gradient\"]\n", - "du_dx = grad_u[\"zonal_gradient\"]\n", - "dv_dy = grad_v[\"meridional_gradient\"]\n", - "\n", - "curl_radial = dv_dx - du_dy\n", - "div_radial = du_dx + dv_dy\n", + "# Compute curl and divergence (default: scaled by sphere_radius)\n", + "curl_radial = u_radial.curl(v_radial)\n", + "div_radial = u_radial.divergence(v_radial)\n", "\n", "print(\n", " f\"Radial field curl range: [{curl_radial.min().values:.2e}, {curl_radial.max().values:.2e}]\"\n", @@ -863,7 +830,7 @@ "Let's verify some fundamental vector calculus identities using our computed fields:\n", "\n", "### Identity 1: Curl of Gradient (Discretization Residual)\n", - "In the continuous setting: ∇ × (∇φ) = 0. UXarray's finite-volume operators are not mimetic, so this holds only approximately. The residual below is discretization error, not a numerical bug. With the default radius scaling it is ~O(10⁻¹³) for this φ ~ O(1) field (the curl stencil applies a 1/radius² factor); on the unit sphere it is ~O(1–10). It shrinks with grid refinement.\n" + "In the continuous setting: ∇ × (∇φ) = 0. UXarray's finite-volume operators are not mimetic, so this holds only approximately. The residual below is discretization error, not a numerical bug. With the default radius scaling it is ~O(10⁻¹⁵) for this φ ~ O(1) field (the curl stencil applies a 1/radius² factor); on the unit sphere it is ~O(0.1). It shrinks with grid refinement.\n" ] }, { From 3b3399ef9941572fe0797788d3eeb355c9d5e352 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Fri, 11 Sep 2026 11:29:17 -0500 Subject: [PATCH 3/6] Fix solid-body vorticity sign, define curl notation, reconcile residual magnitudes --- docs/user-guide/vector_calculus.ipynb | 40 ++++++++++++++++----------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index 3cce25140..710249c52 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/docs/user-guide/vector_calculus.ipynb @@ -228,13 +228,16 @@ "\n", "### Background\n", "\n", - "The curl of a vector field **F** = (u, v) measures the local rotation or circulation. In 2D it produces a scalar field. On a sphere of radius $a$ it carries a metric term that the planar formula leaves out:\n", + "The curl of a vector field **F** = (u, v) measures local rotation. Here $u$, $v$ are the zonal and meridional components, $\\lambda$ is longitude, $\\varphi$ is latitude, $a$ is the sphere radius (`uxgrid.sphere_radius`), and $dx = a\\cos\\varphi\\,d\\lambda$, $dy = a\\,d\\varphi$ are local distances. On a sphere,\n", "\n", "$$\n", - "\\text{curl}(\\mathbf{F}) = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y} + \\frac{u \\tan\\varphi}{a}\n", + "\\text{curl}(\\mathbf{F}) = \\frac{1}{a\\cos\\varphi}\\left[\\frac{\\partial v}{\\partial \\lambda} - \\frac{\\partial (u\\cos\\varphi)}{\\partial \\varphi}\\right]\n", + "= \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y} + \\frac{u \\tan\\varphi}{a}\n", "$$\n", "\n", - "That last term comes from the convergence of the meridians and vanishes only on the equator. Dropping it costs a factor of two on solid-body rotation, so `curl()` includes it. Building the same quantity by hand out of `gradient()` components gives the planar answer instead.\n", + "The product rule on $\\partial(u\\cos\\varphi)/\\partial\\varphi$ produces the last term. It is a metric term: longitude lines converge toward the poles, so a zonal flow rotates even when $\\partial u/\\partial y = 0$. It is zero on the equator. In Example 3 it is as large as the derivative term. See Holton and Hakim, *An Introduction to Dynamic Meteorology*, Chapter 4.\n", + "\n", + "`u.curl(v)` evaluates the full expression; `v.gradient()[\"zonal_gradient\"] - u.gradient()[\"meridional_gradient\"]` gives only the first two terms.\n", "\n", "- **Positive curl**: Counter-clockwise rotation\n", "- **Negative curl**: Clockwise rotation\n", @@ -246,7 +249,8 @@ "\n", "| **Input** | **Usage** | **Output** |\n", "| ---------------------------- | :---------------------: | -------------------- |\n", - "| Vector field (u, v) | `u.curl(v)` | Scalar curl field |\n" + "| Vector field (u, v) | `u.curl(v)` | Scalar curl field |\n", + "" ] }, { @@ -256,9 +260,10 @@ "source": [ "### Constant Fields (Mathematical Validation)\n", "\n", - "On a plane the curl of a constant vector field is zero everywhere. On a sphere it is not: with $u$ constant, both derivative terms drop out and the metric term survives, leaving $\\zeta = u\\tan\\varphi/a$ exactly.\n", + "On a plane the curl of a constant field is zero. On a sphere the derivative terms vanish but the metric term does not, so `u.curl(v)` $= u\\tan\\varphi/a$ exactly.\n", "\n", - "This grid subset spans about ±2° of latitude, so $\\tan\\varphi/a$ stays below $5.5\\times10^{-9}\\ \\mathrm{m^{-1}}$ and the result looks like round-off — but it is a real signal, and we can check it against the closed form.\n" + "This subset spans about ±2° of latitude, so $\\tan\\varphi/a < 5.5\\times10^{-9}\\ \\mathrm{m^{-1}}$ and the result looks like round-off. It is a real signal; the cell below checks it against $u\\tan\\varphi/a$ from the face latitudes.\n", + "" ] }, { @@ -331,7 +336,8 @@ "\n", "In the continuous setting, curl(∇φ) = 0 exactly for any scalar field φ. On an unstructured mesh, however, gradient and curl are independent finite-volume stencils that do not form a discrete de Rham complex, so the identity holds only approximately. The residual is a **discretization error** — not a bug — and shrinks with grid refinement.\n", "\n", - "The *magnitude* of the residual depends on the units. By default `gradient()`/`curl()` scale by `uxgrid.sphere_radius` (Earth ≈ 6.37×10⁶ m), so each derivative picks up a factor of 1/radius. `curl(∇φ)` applies the gradient stencil twice and therefore carries a factor of 1/radius² (≈ 2×10⁻¹⁴). For the Gaussian field here (φ ~ O(1)) the scaled residual is ~O(10⁻¹⁵); on the unit sphere (`scale_by_radius=False`) the same residual is ~O(0.1). Either way it shrinks with refinement.\n" + "The *magnitude* of the residual depends on the units. By default `gradient()`/`curl()` scale by `uxgrid.sphere_radius` (Earth ≈ 6.37×10⁶ m), so each derivative picks up a factor of 1/radius. `curl(∇φ)` applies the gradient stencil twice and therefore carries a factor of 1/radius² (≈ 2.5×10⁻¹⁴). For the Gaussian field here (φ ~ O(1)) the scaled residual is a few ×10⁻¹⁵; on the unit sphere (`scale_by_radius=False`) the same residual is ~O(0.1). Either way it shrinks with refinement.\n", + "" ] }, { @@ -450,7 +456,7 @@ "source": [ "### Example 3: Relative Vorticity from Solid-Body Rotation\n", "\n", - "The 2D curl on the sphere is **relative vorticity** $\\zeta = \\partial v/\\partial x - \\partial u/\\partial y$, a quantity meteorologists and oceanographers care about every day. With the default `scale_by_radius=True`, `u.curl(v)` returns $\\zeta$ in physical units of $s^{-1}$.\n", + "The curl on the sphere is **relative vorticity** $\\zeta$. With the default `scale_by_radius=True`, `u.curl(v)` returns $\\zeta$ in $s^{-1}$.\n", "\n", "A clean analytical check is **solid-body rotation about the polar axis** with angular speed $\\Omega$:\n", "\n", @@ -458,13 +464,11 @@ "u(\\varphi) = \\Omega R \\cos\\varphi, \\qquad v = 0\n", "$$\n", "\n", - "For this flow the relative vorticity on a sphere of radius $R$ is\n", - "\n", "$$\n", - "\\zeta = -\\frac{1}{R\\cos\\varphi}\\frac{\\partial(u\\cos\\varphi)}{\\partial \\varphi} = -2\\Omega \\sin\\varphi.\n", + "\\zeta = -\\frac{1}{R\\cos\\varphi}\\frac{\\partial(\\Omega R\\cos^2\\varphi)}{\\partial \\varphi} = 2\\Omega \\sin\\varphi\n", "$$\n", "\n", - "We construct the field on the existing MPAS subset and compare `u.curl(v)` against this closed form." + "This is the planetary vorticity $f$. The derivative term and the metric term each contribute $\\Omega\\sin\\varphi$, so the planar formula would return half the answer. We compare `u.curl(v)` against $2\\Omega\\sin\\varphi$ on the MPAS subset." ] }, { @@ -497,7 +501,7 @@ ")\n", "\n", "zeta = u_sbr.curl(v_sbr)\n", - "zeta_analytic = -2.0 * OMEGA * np.sin(lat_rad)\n", + "zeta_analytic = 2.0 * OMEGA * np.sin(lat_rad)\n", "\n", "finite = np.isfinite(zeta.values)\n", "err = np.abs(zeta.values[finite] - zeta_analytic[finite])\n", @@ -510,7 +514,9 @@ "print(\n", " f\"analytic zeta range: [{zeta_analytic.min():.3e}, {zeta_analytic.max():.3e}] 1/s\"\n", ")\n", - "print(f\"max |error|: {err.max():.3e} 1/s (≈ {err.max() / (2 * OMEGA):.1%} of 2Omega)\")" + "print(\n", + " f\"max |error|: {err.max():.3e} 1/s (≈ {err.max() / np.abs(zeta_analytic).max():.2%} of max |zeta|)\"\n", + ")" ] }, { @@ -661,7 +667,8 @@ "source": [ "### Example 2: Divergence of Vortex Field\n", "\n", - "Pure rotation is incompressible, so the divergence should be ≈ 0. The metric term is not identically zero here — $v$ varies across the subset — so the residual sits an order of magnitude above the pure discretization error of the planar formula, and is still negligible against the ~10⁻⁵ signal in the radial case below.\n" + "Pure rotation is incompressible, so the divergence should be ≈ 0. The metric term $-v\\tan\\varphi/a$ is not identically zero here, so the residual (~4×10⁻⁹) is an order of magnitude above the derivative terms alone (~3×10⁻¹⁰) and four orders below the ~10⁻⁵ signal in the radial case below.\n", + "" ] }, { @@ -830,7 +837,8 @@ "Let's verify some fundamental vector calculus identities using our computed fields:\n", "\n", "### Identity 1: Curl of Gradient (Discretization Residual)\n", - "In the continuous setting: ∇ × (∇φ) = 0. UXarray's finite-volume operators are not mimetic, so this holds only approximately. The residual below is discretization error, not a numerical bug. With the default radius scaling it is ~O(10⁻¹⁵) for this φ ~ O(1) field (the curl stencil applies a 1/radius² factor); on the unit sphere it is ~O(0.1). It shrinks with grid refinement.\n" + "In the continuous setting: ∇ × (∇φ) = 0. UXarray's finite-volume operators are not mimetic, so this holds only approximately. The residual below is discretization error, not a numerical bug. With the default radius scaling it is a few ×10⁻¹⁵ for this φ ~ O(1) field (the curl stencil applies a 1/radius² factor); on the unit sphere it is ~O(0.1). It shrinks with grid refinement.\n", + "" ] }, { From 88b301629fb721e0e89e3af6c204038c5923cdbe Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 24 Sep 2026 11:32:41 -0400 Subject: [PATCH 4/6] Clarify spherical curl reference and comparison --- docs/user-guide/vector_calculus.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index 710249c52..22456c8fc 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/docs/user-guide/vector_calculus.ipynb @@ -235,7 +235,7 @@ "= \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y} + \\frac{u \\tan\\varphi}{a}\n", "$$\n", "\n", - "The product rule on $\\partial(u\\cos\\varphi)/\\partial\\varphi$ produces the last term. It is a metric term: longitude lines converge toward the poles, so a zonal flow rotates even when $\\partial u/\\partial y = 0$. It is zero on the equator. In Example 3 it is as large as the derivative term. See Holton and Hakim, *An Introduction to Dynamic Meteorology*, Chapter 4.\n", + "The product rule on $\\partial(u\\cos\\varphi)/\\partial\\varphi$ produces the last term. It is a metric term: longitude lines converge toward the poles, so a zonal flow rotates even when $\\partial u/\\partial y = 0$. It is zero on the equator. In Example 3 it is as large as the derivative term. For the spherical-coordinate geometry behind this term, see Holton and Hakim, *An Introduction to Dynamic Meteorology*, Section 2.3, especially Eqs. (2.9) and (2.12).\n", "\n", "`u.curl(v)` evaluates the full expression; `v.gradient()[\"zonal_gradient\"] - u.gradient()[\"meridional_gradient\"]` gives only the first two terms.\n", "\n", @@ -512,10 +512,10 @@ " f\"computed zeta range: [{zeta.values[finite].min():.3e}, {zeta.values[finite].max():.3e}] 1/s\"\n", ")\n", "print(\n", - " f\"analytic zeta range: [{zeta_analytic.min():.3e}, {zeta_analytic.max():.3e}] 1/s\"\n", + " f\"analytic zeta range: [{zeta_analytic[finite].min():.3e}, {zeta_analytic[finite].max():.3e}] 1/s\"\n", ")\n", "print(\n", - " f\"max |error|: {err.max():.3e} 1/s (≈ {err.max() / np.abs(zeta_analytic).max():.2%} of max |zeta|)\"\n", + " f\"max |error|: {err.max():.3e} 1/s (≈ {err.max() / np.abs(zeta_analytic[finite]).max():.2%} of max |zeta|)\"\n", ")" ] }, From 5e149c020710b9038fc1d0d0fd5b3dd3055baa22 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 24 Sep 2026 16:44:26 -0400 Subject: [PATCH 5/6] Ground the solid-body rotation example --- docs/user-guide/vector_calculus.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index 22456c8fc..959a6a70d 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/docs/user-guide/vector_calculus.ipynb @@ -456,9 +456,9 @@ "source": [ "### Example 3: Relative Vorticity from Solid-Body Rotation\n", "\n", - "The curl on the sphere is **relative vorticity** $\\zeta$. With the default `scale_by_radius=True`, `u.curl(v)` returns $\\zeta$ in $s^{-1}$.\n", + "The curl on the sphere is **relative vorticity** $\\zeta$: the local spinning motion of the flow relative to Earth's surface. With the default `scale_by_radius=True`, `u.curl(v)` returns $\\zeta$ in $s^{-1}$.\n", "\n", - "A clean analytical check is **solid-body rotation about the polar axis** with angular speed $\\Omega$:\n", + "A clean analytical check is **solid-body rotation about the polar axis**, in which every point circles the axis with the same angular speed $\\Omega$. Here $R$ is Earth's radius, $\\varphi$ is latitude, and $u$ and $v$ are the eastward and northward speeds. The eastward speed decreases toward the poles because a point there traces a smaller circle:\n", "\n", "$$\n", "u(\\varphi) = \\Omega R \\cos\\varphi, \\qquad v = 0\n", @@ -468,7 +468,7 @@ "\\zeta = -\\frac{1}{R\\cos\\varphi}\\frac{\\partial(\\Omega R\\cos^2\\varphi)}{\\partial \\varphi} = 2\\Omega \\sin\\varphi\n", "$$\n", "\n", - "This is the planetary vorticity $f$. The derivative term and the metric term each contribute $\\Omega\\sin\\varphi$, so the planar formula would return half the answer. We compare `u.curl(v)` against $2\\Omega\\sin\\varphi$ on the MPAS subset." + "The result is the **planetary vorticity** $f$, the local vertical component of the rotation associated with Earth. The derivative term and the metric term each contribute $\\Omega\\sin\\varphi$, so a flat-plane formula, which omits the sphere's curvature, would return half the answer. We compare `u.curl(v)` against $2\\Omega\\sin\\varphi$ on the MPAS subset." ] }, { From d13add8a1c67ff720624673f97cd80cc1e46e77e Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Thu, 24 Sep 2026 17:27:45 -0400 Subject: [PATCH 6/6] Cite full spherical curl equation --- docs/user-guide/vector_calculus.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index 959a6a70d..ebbabcb58 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/docs/user-guide/vector_calculus.ipynb @@ -235,7 +235,7 @@ "= \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y} + \\frac{u \\tan\\varphi}{a}\n", "$$\n", "\n", - "The product rule on $\\partial(u\\cos\\varphi)/\\partial\\varphi$ produces the last term. It is a metric term: longitude lines converge toward the poles, so a zonal flow rotates even when $\\partial u/\\partial y = 0$. It is zero on the equator. In Example 3 it is as large as the derivative term. For the spherical-coordinate geometry behind this term, see Holton and Hakim, *An Introduction to Dynamic Meteorology*, Section 2.3, especially Eqs. (2.9) and (2.12).\n", + "The product rule on $\\partial(u\\cos\\varphi)/\\partial\\varphi$ produces the last term. It is a metric term: longitude lines converge toward the poles, so a zonal flow rotates even when $\\partial u/\\partial y = 0$. It is zero on the equator. In Example 3 it is as large as the derivative term. The [NCAR Command Language relative-vorticity documentation](https://www.ncl.ucar.edu/Document/Functions/Built-in/uv2vr_cfd.shtml) gives this full expanded equation and cites Bluestein, *Synoptic-Dynamic Meteorology in Midlatitudes* (1992), pp. 113–114.\n", "\n", "`u.curl(v)` evaluates the full expression; `v.gradient()[\"zonal_gradient\"] - u.gradient()[\"meridional_gradient\"]` gives only the first two terms.\n", "\n",