Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
151 changes: 63 additions & 88 deletions docs/user-guide/vector_calculus.ipynb
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,19 @@
"\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 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}) = \\nabla \\times \\mathbf{F} = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y}\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",
"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",
"- **Negative curl**: Clockwise rotation\n",
"- **Zero curl**: No local rotation (irrotational flow)\n",
"\n",
"### Usage\n",
Expand All @@ -244,7 +249,8 @@
"\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 +260,10 @@
"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 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 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",
""
]
},
{
Expand All @@ -268,12 +277,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 +288,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 +315,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 +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² (≈ 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.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",
""
]
},
{
Expand All @@ -342,25 +347,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 +406,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 @@ -462,21 +456,19 @@
"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",
"$$\n",
"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."
]
},
{
Expand Down Expand Up @@ -509,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",
Expand All @@ -522,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",
")"
]
},
{
Expand All @@ -546,10 +540,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 +556,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 +566,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 +580,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 +589,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 +615,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 +642,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 +667,8 @@
"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 $-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",
""
]
},
{
Expand All @@ -692,20 +678,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 +721,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 +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(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 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",
""
]
},
{
Expand Down
2 changes: 1 addition & 1 deletion uxarray/core/dataarray.py
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