Vector Calculus

Vector Calculus#

Computing and working with vector calculus on unstructured grids is essential for analyzing scalar and vector fields in geoscience. This notebook demonstrates UXarray’s finite-volume implementations and verifies expected identities.

We will showcase:

Gradient of face-centered scalar fields (zonal and meridional components)
Curl of vector fields (including curl of a gradient ≈ 0 and a synthetic vortex)
Divergence of vector fields (including Laplacian as ∇·∇φ, divergence of a vortex ≈ 0, and radial expansion > 0)
Vector calculus identities verification

import holoviews as hv
import numpy as np

import uxarray as ux

hv.extension("bokeh")

1. Gradient (∇φ)#

Background#

The gradient of a scalar field φ gives both the direction of steepest increase and the rate of change in that direction. On unstructured grids, we use the Green–Gauss theorem:

\[ \int_V \nabla\phi \, dV = \oint_{\partial V} \phi \, dS \]

Implementation#

In a finite-volume context, the gradient of a scalar field \(\phi\) is obtained by summing fluxes across each cell face and dividing by the cell’s volume.

Input	Usage	Output
Scalar field \(\phi\)	`φ.gradient()`	Vector field \(\nabla\phi\)

Finite-volume discretization#

\[ \int_V \nabla\phi \, dV = \oint_{\partial V} \phi \, dS \]

Discrete gradient at cell center \(C^*\)#

\[ \nabla\phi(C^*) \;\approx\; \frac{1}{\mathrm{Vol}(C^*)} \sum_{f\in\partial C^*} \left( \frac{\phi(C_i) + \phi(C_j)}{2} \right) \;l_{ij}\;\mathbf{n}_{ij} \]

Plotting Gradients#

To visualize gradients, we represent them as vector fields and overlay them on the original scalar data.

def plot_gradient_vectors(uxda_grad, **kwargs):
    """
    Plots gradient vectors using HoloViews
    """
    uxgrid = uxda_grad.uxgrid
    mag = np.hypot(uxda_grad.zonal_gradient, uxda_grad.meridional_gradient)
    angle = np.arctan2(uxda_grad.meridional_gradient, uxda_grad.zonal_gradient)

    return hv.VectorField(
        (uxgrid.face_lon, uxgrid.face_lat, angle, mag), **kwargs
    ).opts(magnitude="Magnitude")

# Overlay the gradient vector field on top of the original data variable
p1 = (
    uxds["face_lat"].plot(cmap="Oranges", aspect=1) * plot_gradient_vectors(grad_lat)
).opts(title="∇ Cell Latitudes")
p2 = (
    uxds["face_lon"].plot(cmap="Oranges", aspect=1) * plot_gradient_vectors(grad_lon)
).opts(title="∇ Cell Longitudes")
p3 = (
    uxds["gaussian"].plot(cmap="Oranges", aspect=1) * plot_gradient_vectors(grad_gauss)
).opts(title="∇ Gaussian")
p4 = (
    uxds["inverse_gaussian"].plot(cmap="Oranges", aspect=1)
    * plot_gradient_vectors(grad_inv_gauss)
).opts(title="∇ Inverse Gaussian")

# Compose all four plots in a 2 column layout
(p1 + p2 + p3 + p4).cols(2).opts(shared_axes=False)

2. Curl (∇ × F)#

Background#

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:

\[ \text{curl}(\mathbf{F}) = \nabla \times \mathbf{F} = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \]

Positive curl: Counter-clockwise rotation
Negative curl: Clockwise rotation
Zero curl: No local rotation (irrotational flow)

Usage#

Curl can be computed using the UxDataArray.curl() method on vector field components:

Input	Usage	Output
Vector field (u, v)	`u.curl(v)`	Scalar curl field

Constant Fields (Mathematical Validation)#

The curl of a constant vector field should be zero everywhere (within numerical precision).

# Constant vector field: u=1, v=2 (face-centered)
u_constant = uxds["face_lat"] * 0 + 1.0
v_constant = uxds["face_lat"] * 0 + 2.0

# Compute partials via gradient
grad_u = u_constant.gradient()
grad_v = v_constant.gradient()
du_dy = grad_u["meridional_gradient"]
dv_dx = grad_v["zonal_gradient"]
curl_constant = dv_dx - du_dy

finite = np.isfinite(curl_constant.values)
vals = curl_constant.values[finite]
print(
    f"Total faces: {curl_constant.size}, interior: {vals.size}, boundary NaNs: {np.isnan(curl_constant.values).sum()}"
)
if vals.size:
    print(f"Finite curl range: [{vals.min():.2e}, {vals.max():.2e}]")
    print(
        f"Max |curl|: {np.abs(vals).max():.2e}, Mean |curl|: {np.abs(vals).mean():.2e}"
    )

Total faces: 195, interior: 147, boundary NaNs: 48
Finite curl range: [0.00e+00, 0.00e+00]
Max |curl|: 0.00e+00, Mean |curl|: 0.00e+00

Gaussian Fields#

# Use Gaussian fields as vector components
u_gauss = uxds["gaussian"]
v_gauss = uxds["inverse_gaussian"]

# Compute partials via gradient
grad_u = u_gauss.gradient()
grad_v = v_gauss.gradient()
du_dy = grad_u["meridional_gradient"]
dv_dx = grad_v["zonal_gradient"]
curl_gauss = dv_dx - du_dy

finite = np.isfinite(curl_gauss.values)
vals = curl_gauss.values[finite]
print(
    f"Total faces: {curl_gauss.size}, interior: {vals.size}, boundary NaNs: {np.isnan(curl_gauss.values).sum()}"
)
if vals.size:
    print(f"Finite curl range: [{vals.min():.6f}, {vals.max():.6f}]")
    print(f"Mean curl (finite): {vals.mean():.6f}")

Total faces: 195, interior: 147, boundary NaNs: 48
Finite curl range: [-142.311439, 142.830187]
Mean curl (finite): -0.320411

Example 1: Curl of Gradient Field#

A fundamental vector identity states that the curl of any gradient field should be zero (conservative field). Let’s verify this:

# Extract gradient components
u_component = grad_gauss.zonal_gradient
v_component = grad_gauss.meridional_gradient

# Compute partial derivatives via gradient()
grad_u = u_component.gradient()
grad_v = v_component.gradient()
du_dy = grad_u["meridional_gradient"]
dv_dx = grad_v["zonal_gradient"]

# Curl = ∂v/∂x - ∂u/∂y
curl_of_gradient = dv_dx - du_dy

print(
    f"Curl of gradient range: [{curl_of_gradient.min().values:.2e}, {curl_of_gradient.max().values:.2e}]"
)
print(f"Mean absolute curl: {abs(curl_of_gradient).mean().values:.2e}")

curl_plot = curl_of_gradient.plot(cmap="RdBu_r", aspect=1, clim=(-1e-10, 1e-10)).opts(
    title="Curl of Gradient Field (Should ≈ 0)", colorbar=True
)
curl_plot

Curl of gradient range: [-4.51e+00, 4.35e+00]
Mean absolute curl: 1.76e+00

Example 2: Synthetic Vortex Field#

To better demonstrate curl, let’s create a synthetic rotating vector field (vortex):

# Get face coordinates
face_lon = uxds.uxgrid.face_lon.values
face_lat = uxds.uxgrid.face_lat.values

# Center the coordinates
lon_center = np.mean(face_lon)
lat_center = np.mean(face_lat)

x = face_lon - lon_center
y = face_lat - lat_center

# Create a vortex: u = -y, v = x (pure rotation)
u_vortex_data = -y
v_vortex_data = x

# Create UxDataArrays
u_vortex = ux.UxDataArray(
    u_vortex_data, dims=["n_face"], uxgrid=uxds.uxgrid, name="u_vortex"
)
v_vortex = ux.UxDataArray(
    v_vortex_data, dims=["n_face"], uxgrid=uxds.uxgrid, name="v_vortex"
)

# Compute curl via gradients
grad_u = u_vortex.gradient()
grad_v = v_vortex.gradient()
du_dy = grad_u["meridional_gradient"]
dv_dx = grad_v["zonal_gradient"]
curl_vortex = dv_dx - du_dy

print(
    f"Vortex curl range: [{curl_vortex.min().values:.2f}, {curl_vortex.max().values:.2f}]"
)
print("Expected curl for pure rotation: ~2.0")

Vortex curl range: [343.76, 343.85]
Expected curl for pure rotation: ~2.0

# Plot the vortex vector field and its curl
def plot_vector_field(u, v, **kwargs):
    """Plot vector field using HoloViews"""
    uxgrid = u.uxgrid
    mag = np.hypot(u, v)
    angle = np.arctan2(v, u)

    return hv.VectorField(
        (uxgrid.face_lon, uxgrid.face_lat, angle, mag), **kwargs
    ).opts(magnitude="Magnitude")


# Vector field plot
vortex_vectors = plot_vector_field(u_vortex, v_vortex).opts(
    title="Synthetic Vortex Field", aspect=1, arrow_heads=True
)

# Curl magnitude plot
curl_magnitude = curl_vortex.plot(cmap="Reds", aspect=1).opts(
    title="Curl of Vortex Field", colorbar=True
)

(vortex_vectors + curl_magnitude).opts(shared_axes=False)

3. Divergence (∇ · F)#

Background#

The divergence of a vector field F = (u, v) measures the local expansion or contraction of the field:

\[ \text{div}(\mathbf{F}) = \nabla \cdot \mathbf{F} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \]

Positive divergence: Expansion (source)
Negative divergence: Contraction (sink)
Zero divergence: Incompressible flow

Usage#

Divergence can be computed using the UxDataArray.divergence() method:

Input	Usage	Output
Vector field (u, v)	`u.divergence(v)`	Scalar divergence field

Constant Fields (Mathematical Validation)#

The divergence of a constant vector field should be zero everywhere (within numerical precision).

# Constant vector field: u=1, v=2 (face-centered)
u_constant = uxds["face_lat"] * 0 + 1.0
v_constant = uxds["face_lat"] * 0 + 2.0

# Compute partials via gradient
grad_u = u_constant.gradient()
grad_v = v_constant.gradient()
du_dx = grad_u["zonal_gradient"]
dv_dy = grad_v["meridional_gradient"]
div_constant = du_dx + dv_dy

finite = np.isfinite(div_constant.values)
vals = div_constant.values[finite]
print(
    f"Total faces: {div_constant.size}, interior: {vals.size}, boundary NaNs: {np.isnan(div_constant.values).sum()}"
)
if vals.size:
    print(f"Finite divergence range: [{vals.min():.2e}, {vals.max():.2e}]")
    print(f"Max |div|: {np.abs(vals).max():.2e}, Mean |div|: {np.abs(vals).mean():.2e}")

Total faces: 195, interior: 147, boundary NaNs: 48
Finite divergence range: [0.00e+00, 0.00e+00]
Max |div|: 0.00e+00, Mean |div|: 0.00e+00

Gaussian Fields#

# Use Gaussian fields as vector components
u_gauss = uxds["gaussian"]
v_gauss = uxds["inverse_gaussian"]

# Compute partials via gradient
grad_u = u_gauss.gradient()
grad_v = v_gauss.gradient()
du_dx = grad_u["zonal_gradient"]
dv_dy = grad_v["meridional_gradient"]
div_gauss = du_dx + dv_dy

finite = np.isfinite(div_gauss.values)
vals = div_gauss.values[finite]
print(
    f"Total faces: {div_gauss.size}, interior: {vals.size}, boundary NaNs: {np.isnan(div_gauss.values).sum()}"
)
if vals.size:
    print(f"Finite divergence range: [{vals.min():.6f}, {vals.max():.6f}]")
    print(f"Mean divergence (finite): {vals.mean():.6f}")

Total faces: 195, interior: 147, boundary NaNs: 48
Finite divergence range: [-141.595927, 143.313236]
Mean divergence (finite): 0.974066

Example 1: Divergence of Gradient Field#

# Compute divergence of the gradient (Laplacian) via partials
# ∇²φ = ∂(∇φ_x)/∂x + ∂(∇φ_y)/∂y
grad_gauss_u = grad_gauss["zonal_gradient"]
grad_gauss_v = grad_gauss["meridional_gradient"]

gxu = grad_gauss_u.gradient()["zonal_gradient"]  # ∂u/∂x
gyv = grad_gauss_v.gradient()["meridional_gradient"]  # ∂v/∂y

div_of_gradient = gxu + gyv

print(
    f"Divergence of gradient range: [{div_of_gradient.min().values:.2e}, {div_of_gradient.max().values:.2e}]"
)
print("This is the Laplacian (∇²) of the original gaussian field")

div_plot = div_of_gradient.plot(cmap="RdBu_r", aspect=1).opts(
    title="Divergence of Gradient (Laplacian)", colorbar=True
)
div_plot

Divergence of gradient range: [-5.45e+04, 1.03e+03]
This is the Laplacian (∇²) of the original gaussian field

Example 2: Divergence of Vortex Field#

Pure rotation should have zero divergence (incompressible):

# Compute divergence of the vortex via gradients: div = ∂u/∂x + ∂v/∂y
grad_u = u_vortex.gradient()
grad_v = v_vortex.gradient()
du_dx = grad_u["zonal_gradient"]
dv_dy = grad_v["meridional_gradient"]
div_vortex = du_dx + dv_dy

print(
    f"Vortex divergence range: [{div_vortex.min().values:.2e}, {div_vortex.max().values:.2e}]"
)
print(f"Mean absolute divergence: {abs(div_vortex).mean().values:.2e}")
print("Pure rotation should have zero divergence")

div_vortex_plot = div_vortex.plot(cmap="RdBu_r", aspect=1, clim=(-1e-10, 1e-10)).opts(
    title="Divergence of Vortex (Should ≈ 0)", colorbar=True
)
div_vortex_plot

Vortex divergence range: [-5.41e-03, 4.81e-03]
Mean absolute divergence: 1.72e-03
Pure rotation should have zero divergence

Example 3: Radial Expansion Field#

Let’s create a field that expands radially outward to demonstrate positive divergence:

# Create radial expansion field: u = x, v = y
u_radial_data = x
v_radial_data = y

u_radial = ux.UxDataArray(
    u_radial_data, dims=["n_face"], uxgrid=uxds.uxgrid, name="u_radial"
)
v_radial = ux.UxDataArray(
    v_radial_data, dims=["n_face"], uxgrid=uxds.uxgrid, name="v_radial"
)

# Compute curl and divergence via gradients
grad_u = u_radial.gradient()
grad_v = v_radial.gradient()
du_dy = grad_u["meridional_gradient"]
dv_dx = grad_v["zonal_gradient"]
du_dx = grad_u["zonal_gradient"]
dv_dy = grad_v["meridional_gradient"]

curl_radial = dv_dx - du_dy
div_radial = du_dx + dv_dy

print(
    f"Radial field curl range: [{curl_radial.min().values:.2e}, {curl_radial.max().values:.2e}]"
)
print(
    f"Radial field divergence range: [{div_radial.min().values:.2f}, {div_radial.max().values:.2f}]"
)
print("Expected: curl ≈ 0, divergence ≈ 2 (pure expansion)")

Radial field curl range: [-4.81e-03, 5.41e-03]
Radial field divergence range: [343.76, 343.85]
Expected: curl ≈ 0, divergence ≈ 2 (pure expansion)

# Plot radial expansion field and its divergence
radial_vectors = plot_vector_field(u_radial, v_radial).opts(
    title="Radial Expansion Field", aspect=1, arrow_heads=True
)

radial_div = div_radial.plot(cmap="Reds", aspect=1).opts(
    title="Divergence of Radial Field", colorbar=True
)

(radial_vectors + radial_div).opts(shared_axes=False)

4. Vector Calculus Identities#

Let’s verify some fundamental vector calculus identities using our computed fields:

Identity 1: Curl of Gradient is Zero#

For any scalar field φ: ∇ × (∇φ) = 0

# We already computed this above
max_curl_grad = np.abs(curl_of_gradient).max().values
print(f"Maximum |curl(∇φ)|: {max_curl_grad:.2e}")
print(f"Identity verified: {max_curl_grad < 1e-10}")

Maximum |curl(∇φ)|: 4.51e+00
Identity verified: False

Identity 2: Divergence of Curl is Zero#

For any vector field F: ∇ · (∇ × F) = 0

Note: This identity applies to 3D vector fields. In 2D, curl produces a scalar, so we can’t directly compute its divergence.

Identity 3: Properties of Special Fields#

# Summary of field properties
print("Field Properties Summary:")
print("=" * 50)
print("Gradient field:")
print(f"  - Curl: {np.abs(curl_of_gradient).max().values:.2e} (≈ 0, conservative)")
print(
    f"  - Divergence: {div_of_gradient.min().values:.2e} to {div_of_gradient.max().values:.2e} (Laplacian)"
)
print()
print("Vortex field (pure rotation):")
print(f"  - Curl: {curl_vortex.mean().values:.2f} (≈ 2, rotational)")
print(f"  - Divergence: {np.abs(div_vortex).max().values:.2e} (≈ 0, incompressible)")
print()
print("Radial field (pure expansion):")
print(f"  - Curl: {np.abs(curl_radial).max().values:.2e} (≈ 0, irrotational)")
print(f"  - Divergence: {div_radial.mean().values:.2f} (≈ 2, expanding)")

Field Properties Summary:
==================================================
Gradient field:
  - Curl: 4.51e+00 (≈ 0, conservative)
  - Divergence: -5.45e+04 to 1.03e+03 (Laplacian)

Vortex field (pure rotation):
  - Curl: 343.79 (≈ 2, rotational)
  - Divergence: 5.41e-03 (≈ 0, incompressible)

Radial field (pure expansion):
  - Curl: 5.41e-03 (≈ 0, irrotational)
  - Divergence: 343.79 (≈ 2, expanding)

Vector Calculus

Contents

Vector Calculus#

Data#

1. Gradient (∇φ)#

Background#

Implementation#

Finite-volume discretization#

Discrete gradient at cell center \(C^*\)#

Usage#

Plotting Gradients#

2. Curl (∇ × F)#

Background#

Usage#

Constant Fields (Mathematical Validation)#

Gaussian Fields#

Example 1: Curl of Gradient Field#

Example 2: Synthetic Vortex Field#

3. Divergence (∇ · F)#

Background#

Usage#

Constant Fields (Mathematical Validation)#

Gaussian Fields#

Example 1: Divergence of Gradient Field#

Example 2: Divergence of Vortex Field#

Example 3: Radial Expansion Field#

4. Vector Calculus Identities#

Identity 1: Curl of Gradient is Zero#

Identity 2: Divergence of Curl is Zero#

Identity 3: Properties of Special Fields#