Skip to content

Vector calculus follow-ups from #1663: eager divergence and a user guide fix #1732

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rajeeja wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Vector calculus follow-ups from #1663: eager divergence and a user guide fix #1732

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
609f80e
Keep divergence lazy on dask-backed input
rajeeja Sep 3, 2026
248bb8c
Teach the spherical operators in the vector calculus user guide
rajeeja Sep 3, 2026
4a8be4d
Merge branch 'main' into rajeeja/vector-calculus-followup
rajeeja Sep 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
129 changes: 48 additions & 81 deletions docs/user-guide/vector_calculus.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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`."
]
},
{
Expand Down Expand Up @@ -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",

@Sevans711 Sevans711 Sep 9, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to define the meaning of \varphi. Is it just one of the spherical coordinates?

"$$\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",

@Sevans711 Sevans711 Sep 9, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is very hard for me to parse, please include more details. Maybe it would also help to link to a paper or website that presents/derives the formula, too?

Notes:

  • "convergence of the meridians" - what does this mean?
  • "comes from the convergence of the meridians" - how?
  • "Dropping it costs a factor of two on solid-body rotation" - why?
  • "so curl() includes it" - this could be removed; it was a bit confusing to read. My default assumption is that the uxarray implementation would match whatever formula is presented here.
  • "Building the same quantity by hand out of gradient() components" - what does this mean?
  • "gives the planar answer instead" - what does "the planar answer" mean?

"\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",
Expand All @@ -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"
]
},
{
Expand All @@ -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",

@Sevans711 Sevans711 Sep 9, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

zeta is not defined yet; it would be easier to read this section if not defining a new variable, something like this would work great: "leaving u.curl(v) = … exactly"

"\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"

@Sevans711 Sevans711 Sep 9, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is "the closed form"?

]
},
{
Expand All @@ -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",
Expand All @@ -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}\")"
]
},
{
Expand All @@ -306,21 +310,16 @@
"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",
"print(\n",
" 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}\")"
]
},
{
Expand All @@ -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"

@Sevans711 Sevans711 Sep 9, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why did these numbers change? Code changes in this PR should not affect numerical results from these methods, right?

]
},
{
Expand All @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand All @@ -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"
]
},
{
Expand All @@ -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"
]
},
{
Expand All @@ -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",
Expand All @@ -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}\")"
]
},
{
Expand All @@ -622,21 +609,16 @@
"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",
"print(\n",
" 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}\")"
]
},
{
Expand All @@ -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",
Expand All @@ -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"
]
},
{
Expand All @@ -692,20 +671,15 @@
"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",
")\n",
"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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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"
]
},
{
Expand Down
2 changes: 1 addition & 1 deletion uxarray/core/dataarray.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading