From 5bf1eed657a122511bb72511ba2427f2de439a08 Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Mon, 2 Nov 2026 18:25:35 -0500 Subject: [PATCH 1/8] add UxDataArray.laplacian() --- uxarray/core/dataarray.py | 67 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index a10f4ab1c..75fb2329c 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -1946,6 +1946,73 @@ def divergence( result = result.assign_attrs(attrs) return result + def laplacian(self, *, scale_by_radius: bool = True) -> "UxDataArray": + """ + Computes the Laplacian of this UxDataArray in the horizontal direction + (i.e., along the unstructured grid dimension). + + The formula is: laplacian(q) = div(grad(q)) = ∂²q/∂x² + ∂²q/∂y² - (∂q/∂y)·tan(φ)/a, + where q is a scalar field, φ=latitude, and a=sphere radius. + where x corresponds to the zonal direction and y the meridional direction. + + Parameters + ---------- + scale_by_radius : bool, default=True + Whether to divide unit-sphere derivatives by ``uxgrid.sphere_radius``. + When ``True`` (and the grid's ``sphere_radius`` is set correctly), + the result carries physical units: [q units]/meter^2 + (e.g. result has units of 1/s^2 if q represents velocity). + When ``False`` (or if the grid's ``sphere_radius`` is set to 1.0), + the result caries per-radian^2 units: [q units]/radian^2. + (e.g. result has units of (m/s)/rad^2 if q represents velocity). + + CAUTION: this scaling does not account for any radius increases + due to vertical levels; it simply divides by ``uxgrid.sphere_radius``. + + Returns + ------- + laplacian : UxDataArray + The Laplacian of the scalar field, as a UxDataArray. + Units determined by ``scale_by_radius``. + + Example + ------- + >>> uxds["var"].laplacian() + """ + # Compute the gradient of the scalar field + grad_q_zonal, grad_q_meridional = _compute_gradient( + self, scale_by_radius=scale_by_radius + ) + + # Compute the divergence of the gradient to get the Laplacian + result = grad_q_zonal.divergence( + grad_q_meridional, scale_by_radius=scale_by_radius + ) + + # bookkeeping for units + units = result.attrs["units"] + assert units.endswith("/m") or units.endswith( + "/rad" + ) # else, made a coding error. + u_numer, u_denom = units.rsplit("/", 1) + units = f"{u_numer}/{u_denom}^2" # e.g., input has "K", output has "K/m^2". + + # bookkeeping for name & attrs + result = result.rename("laplacian") + self_longname = self.name if self.name is not None else "unnamed_array" + attrs = { + "long_name": f"Laplacian of {self_longname}", + "units": units, + "derivative": "laplacian", + "description": ( + "Laplacian of scalar field q computed as div(grad(q)) " + "== ∂²q/∂x² + ∂²q/∂y² - (∂q/∂y)·tan(φ)/a, " + "with x zonal; y meridional; φ=latitude; a=sphere radius." + ), + } + result = result.assign_attrs(attrs) + return result + def scalardotgradient(self, v: "UxDataArray", q: "UxDataArray") -> "UxDataArray": """ Compute the dot product between a vector field and the gradient of a scalar field. From 9263dad108c27a8ccd43b1808c8e98a8a6ef9d74 Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Mon, 2 Nov 2026 18:36:57 -0500 Subject: [PATCH 2/8] simple tests for laplacian() (still need to actually test for numerical correctness.) --- test/core/test_vector_calculus.py | 87 +++++++++++++++++++++++++++++++ uxarray/core/dataarray.py | 5 +- 2 files changed, 88 insertions(+), 4 deletions(-) diff --git a/test/core/test_vector_calculus.py b/test/core/test_vector_calculus.py index 1df847216..d80e193a5 100644 --- a/test/core/test_vector_calculus.py +++ b/test/core/test_vector_calculus.py @@ -1012,6 +1012,66 @@ def test_curl_units_and_attributes(self, gridpath, datasetpath): assert curl_field.uxgrid == u_component.uxgrid +class TestLaplacianQuadHex: + + def test_laplacian_output_format(self): + """Tests the output format of laplacian functionality""" + uxds = ux.tutorial.open_dataset("quad-hexagon") + q = uxds['t2m'] + laplacian_da = q.laplacian() + + assert isinstance(laplacian_da, ux.UxDataArray) + assert laplacian_da.name == "laplacian" + assert laplacian_da.attrs["derivative"] == "laplacian" + assert laplacian_da.attrs["long_name"] == "Laplacian of t2m" + assert laplacian_da.attrs["description"].startswith("Laplacian of scalar field") + assert q.sizes == laplacian_da.sizes + + # units attr checks + assert q.attrs["units"] == "K" + assert laplacian_da.attrs["units"] == "K/m^2" + laplacian_da_no_scaling = q.laplacian(scale_by_radius=False) + assert laplacian_da_no_scaling.attrs["units"] == "K/rad^2" + q_no_units = q.copy() + q_no_units.attrs.pop("units") + laplacian_da_no_units_yes_scaling = q_no_units.laplacian() + assert laplacian_da_no_units_yes_scaling.attrs["units"] == "1/m^2" + laplacian_da_no_units_no_scaling = q_no_units.laplacian(scale_by_radius=False) + assert laplacian_da_no_units_no_scaling.attrs["units"] == "1/rad^2" + q_ratio_units = q.assign_attrs({"units": "K/s"}) + laplacian_da_ratio_units = q_ratio_units.laplacian() + assert laplacian_da_ratio_units.attrs["units"] == "(K/s)/m^2" + # blank units input should become 1/m^2 or 1/rad^2, not "/m^2" or "/rad^2" + for blank in ["", " ", " \t ", "\n"]: + q_blank_units = q.assign_attrs({"units": blank}) + laplacian_da_blank_units = q_blank_units.laplacian() + assert laplacian_da_blank_units.attrs["units"] == "1/m^2" + laplacian_da_blank_units_no_scaling = q_blank_units.laplacian(scale_by_radius=False) + assert laplacian_da_blank_units_no_scaling.attrs["units"] == "1/rad^2" + # units of None, False, or 0 should become None/m^2, False/m^2, or 0/m^2, not 1/m^2. + # (similarly for rad^2, but just spot-checking m^2 for brevity.) + for zero in [None, False, 0]: + q_zero_units = q.assign_attrs({"units": zero}) + laplacian_da_zero_units = q_zero_units.laplacian() + assert laplacian_da_zero_units.attrs["units"] == f"{zero}/m^2" + + # other attrs checks - spot check for handling None names: + # use "unnamed_array" in long_name; keep name as "laplacian" + q_noname = q.rename(None) + laplacian_da_noname = q_noname.laplacian() + assert laplacian_da_noname.name == "laplacian" + assert laplacian_da_noname.attrs["long_name"] == "Laplacian of unnamed_array" + + # other attrs checks - spot check that value-describing attrs get removed + # (gradient output format test has more comprehensive checks): + q_with_more_attrs = q.copy() + q_with_more_attrs.attrs["actual_range"] = (0.0, 100.0) + q_with_more_attrs.attrs["some other attr"] = "test" + laplacian_with_more_attrs = q_with_more_attrs.laplacian() + assert "actual_range" not in laplacian_with_more_attrs.attrs + assert "some other attr" in laplacian_with_more_attrs.attrs + + class TestSphericalManufacturedSolutions: """Amplitude checks against closed-form answers on the unit sphere. @@ -1264,6 +1324,33 @@ def test_divergence_more_than_1d(self): resultj = arrj.divergence(arr0j) assert np.allclose(resultj.imag, resultA, atol=0, rtol=1e-9) + def test_laplacian_more_than_1d(self): + ds0 = ux.tutorial.open_dataset('outCSne30-vortex') + arr0 = ds0['psi'] + arr0.uxgrid._ds.attrs['sphere_radius'] = 1.0 + # ^ set 'sphere_radius' to avoid the warning which is irrelevant here: + # "scale_by_radius=True but the grid has no 'sphere_radius' attribute" + times = xr.DataArray([1,3,5,7], coords={'time': [10,30,50,70]}) + levels = xr.DataArray([2,4], coords={'level': [2,4]}) + arr = arr0 * times + levels + assert arr.sizes == {'time': 4, 'level': 2, 'n_face': arr0.sizes['n_face']} + result = arr.laplacian() + assert result.sizes == arr.sizes + # ensure full result matches 1D results at each slice + for it in range(arr.sizes['time']): + for il in range(arr.sizes['level']): + arr_it_il_laplacian = arr.isel(time=it, level=il).laplacian() + assert result.isel(time=it, level=il).equals(arr_it_il_laplacian) + # ensure chunked input stays chunked and produces correct results + arrc = arr.chunk({'time': 2, 'level': -1}) + resultc = arrc.laplacian() + assert resultc.chunks == arrc.chunks + assert resultc.compute().equals(result) + # ensure complex values are allowed :) + arrj = arr * 1j + resultj = arrj.laplacian() + assert np.allclose(resultj.imag, result, atol=0, rtol=1e-9) + def test_scalar_dot_gradient_more_than_1d(self): ds0 = ux.tutorial.open_dataset('outCSne30-vortex') arr0 = ds0['psi'] diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index 75fb2329c..891a5c951 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -1990,10 +1990,7 @@ def laplacian(self, *, scale_by_radius: bool = True) -> "UxDataArray": ) # bookkeeping for units - units = result.attrs["units"] - assert units.endswith("/m") or units.endswith( - "/rad" - ) # else, made a coding error. + units = grad_q_zonal.attrs["units"] # ends with "/m" or "/rad". u_numer, u_denom = units.rsplit("/", 1) units = f"{u_numer}/{u_denom}^2" # e.g., input has "K", output has "K/m^2". From d756bc940ac365b14594060c247b08ce0cbb0507 Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Mon, 2 Nov 2026 18:45:51 -0500 Subject: [PATCH 3/8] fix user-guide docs missing laplacian metric term --- docs/user-guide/vector_calculus.ipynb | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index 78b74f1a5..16fb4a718 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/docs/user-guide/vector_calculus.ipynb @@ -645,7 +645,8 @@ "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", "grad_gauss_u = grad_gauss[\"zonal_gradient\"]\n", "grad_gauss_v = grad_gauss[\"meridional_gradient\"]\n", "\n", From c8ac579ed13b0adb705b33fd51ab2eba7c261ad1 Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Mon, 2 Nov 2026 18:56:31 -0500 Subject: [PATCH 4/8] docs: add UxDataArray.laplacian in api.rst --- docs/api.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/api.rst b/docs/api.rst index 88231ed05..58ed3cd4c 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -529,6 +529,7 @@ Mathematical Operators UxDataArray.divergence UxDataArray.gradient UxDataArray.integrate + UxDataArray.laplacian UxDataArray.scalardotgradient From 2a28f4c5809f118858c7f83c6443fd64bfd29fc8 Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Tue, 3 Nov 2026 16:08:01 -0500 Subject: [PATCH 5/8] numerical tests for laplacian --- test/core/test_vector_calculus.py | 159 ++++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) diff --git a/test/core/test_vector_calculus.py b/test/core/test_vector_calculus.py index d80e193a5..9d9e05918 100644 --- a/test/core/test_vector_calculus.py +++ b/test/core/test_vector_calculus.py @@ -1072,6 +1072,165 @@ def test_laplacian_output_format(self): assert "some other attr" in laplacian_with_more_attrs.attrs +class TestLaplacianDyamondSubset: + + def test_laplacian_constant_field(self, gridpath, datasetpath): + """Laplacian of a constant scalar field reduces to 0. + The metric term from div in div(grad(q)) sees v = ∂q/∂y = 0, so + the result is identically 0 (or NaN if missing boundaries) everywhere. + """ + uxds = ux.open_dataset( + gridpath("mpas", "dyamond-30km", "gradient_grid_subset.nc"), + datasetpath("mpas", "dyamond-30km", "gradient_data_subset.nc") + ) + # some constant fields: + q0 = uxds['face_lat'] * 0 + q1 = q0 + 1.0 + q3 = q0 - 3.0 + + lap0 = q0.laplacian() + lap1 = q1.laplacian() + lap3 = q3.laplacian() + boundaries = lap0.isnull() + assert np.all(boundaries == lap1.isnull()) + assert np.all(boundaries == lap3.isnull()) + + assert np.all(lap0[~boundaries] == 0) + assert np.all(lap1[~boundaries] == 0) + assert np.all(lap3[~boundaries] == 0) + + # changing scale_by_radius (nontrivial because the radius is not 1. + assert uxds.uxgrid.sphere_radius > 1.0 + lap0r = q0.laplacian(scale_by_radius=False) + lap1r = q1.laplacian(scale_by_radius=False) + lap3r = q3.laplacian(scale_by_radius=False) + assert np.all(lap0r[~boundaries] == 0) + assert np.all(lap1r[~boundaries] == 0) + assert np.all(lap3r[~boundaries] == 0) + + +class TestLaplacianSphericalHarmonics: + """Numerical correctness checks for the Laplacian, using spherical harmonics: + laplacian(q) == -l(l+1)/r^2 * q + where q is any linear combination of any spherical harmonics of degree l. + + Linear combinations of harmonics of degree 1 include: x, y, z + Linear combinations of harmonics of degree 1 include: xy, xz, yz, x^2-y^2, 3z^2-r^2 + """ + def grid_and_xyz(self, zoom): + """returns (healpix grid, x, y, z) for this zoom. + x, y, z are from the grid, converted to UxDataArrays. + """ + grid = ux.Grid.from_healpix(zoom=zoom) + grid.sphere_radius = 1.0 # avoid any warnings later about undefined sphere_radius. + zero = ux.UxDataArray(xr.DataArray(np.zeros(grid.n_face), dims="n_face"), uxgrid=grid) + x = grid.face_x + zero + y = grid.face_y + zero + z = grid.face_z + zero + return grid, x, y, z + + def error_stats(self, q, l=1, *, max_lat=70): + """returns (mean, max) of |lap(q) - expected| / max(q), + with expected = -l(l+1)/r^2 * q, and the difference is only computed where |lat| < max_lat. + Restricting to mid-latitudes large tan(lat) near poles. + + mean gives a decent measurement of overall error, while max tells worst-case error. + """ + max_q = q.max().item() + ylm_factor = -l * (l+1) / q.uxgrid.sphere_radius**2 + diff = q.laplacian() - ylm_factor * q + far_from_poles = np.abs(q.uxgrid.face_lat) < max_lat + diff_far_from_poles = diff[far_from_poles] + err_check = np.abs(diff_far_from_poles) + return err_check.mean().item() / max_q, err_check.max().item() / max_q + + def test_laplacian_harmonics_l1_zoom4(self): + """checks absolute errors some spherical harmonics on healpix grid at zoom=4. + The error bounds here were derived empirically; they are not tiny because uxarray's + Laplacian uses a low-order stencil (see test_laplacian_harmonics_l1_convergence). + """ + MEAN_TOL = 0.02 + MAX_TOL = 0.15 + grid, x, y, z = self.grid_and_xyz(zoom=4) + # various combinations of spherical harmonics + fields = [x, y, z, 7*x, y+z, 3*y-10*z, -x-y-z, 2*x-3*y+4*z] + for i, q in enumerate(fields): + mean_err, max_err = self.error_stats(q, l=1) + assert mean_err < MEAN_TOL, f"failure on field i={i}" + assert max_err < MAX_TOL, f"failure on field i={i}" + + def test_laplacian_harmonics_l2_zoom4(self): + """checks absolute errors some spherical harmonics on healpix grid at zoom=4. + The error bounds here were derived empirically; they are not tiny because uxarray's + Laplacian uses a low-order stencil (see test_laplacian_harmonics_l2_convergence). + """ + MEAN_TOL = 0.07 + MAX_TOL = 0.50 + grid, x, y, z = self.grid_and_xyz(zoom=4) + # various combinations of spherical harmonics of degree 2 + fields = [x*y, x*z, y*z, x**2-y**2, 3*z**2-1, 7*x*y, -3*x*z+5*y*z+7*(x**2-y**2)] + for i, q in enumerate(fields): + mean_err, max_err = self.error_stats(q, l=2) + assert mean_err < MEAN_TOL, f"failure on field i={i}" + assert max_err < MAX_TOL, f"failure on field i={i}" + + def test_laplacian_harmonics_l1_convergence(self): + """checks that the Laplacian mean errors shrink as the grid is refined, + for spherical harmonics of degree 1. (Max errors not tested because the low-order stencil + happens to not always cause max errors to decrease as zoom increases. Checking by-eye + reveals these max errors occurs near some of few special vertices where only 3 faces meet.) + """ + zooms = [4, 5, 6, 7] + # expect at least this much improvement in mean error when zoom increases by 1. + MEAN_CONV_TOL = 3.0 + zoom_to_errors = {} # {zoom: list of errors for each field} + for zoom in zooms: + errors = [] + grid, x, y, z = self.grid_and_xyz(zoom=zoom) + fields = [x, y, z, 7*x, y+z] + for q in fields: + mean_err, max_err = self.error_stats(q, l=1) + errors.append(mean_err) + zoom_to_errors[zoom] = errors + for j in range(len(fields)): + for i in range(len(zooms)-1): + zoom1, zoom2 = zooms[i], zooms[i+1] + err1 = zoom_to_errors[zoom1][j] + err2 = zoom_to_errors[zoom2][j] + assert err1 / err2 >= MEAN_CONV_TOL, ( + f"mean error did not shrink enough from zoom {zoom1} to {zoom2} " + f"for field {j}: {err1} -> {err2}" + ) + + def test_laplacian_harmonics_l2_convergence(self): + """checks that the Laplacian mean errors shrink as the grid is refined, + for spherical harmonics of degree 2. (Max errors not tested because the low-order stencil + happens to not always cause max errors to decrease as zoom increases. Checking by-eye + reveals these max errors occurs near some of few special vertices where only 3 faces meet.) + """ + zooms = [4, 5, 6, 7] + # expect at least this much improvement in mean error when zoom increases by 1. + MEAN_CONV_TOL = 2.0 + zoom_to_errors = {} # {zoom: list of errors for each field} + for zoom in zooms: + errors = [] + grid, x, y, z = self.grid_and_xyz(zoom=zoom) + fields = [x*y, x*z, y*z, x**2-y**2, 3*z**2-1] + for q in fields: + mean_err, max_err = self.error_stats(q, l=2) + errors.append(mean_err) + zoom_to_errors[zoom] = errors + for j in range(len(fields)): + for i in range(len(zooms)-1): + zoom1, zoom2 = zooms[i], zooms[i+1] + err1 = zoom_to_errors[zoom1][j] + err2 = zoom_to_errors[zoom2][j] + assert err1 / err2 >= MEAN_CONV_TOL, ( + f"mean error did not shrink enough from zoom {zoom1} to {zoom2} " + f"for field {j}: {err1} -> {err2}" + ) + + class TestSphericalManufacturedSolutions: """Amplitude checks against closed-form answers on the unit sphere. From 960b015e3aeca433a7fadf8dd6318b09bd27cb3d Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Tue, 3 Nov 2026 16:36:05 -0500 Subject: [PATCH 6/8] add laplacian examples to vector_calculus.ipynb --- docs/user-guide/vector_calculus.ipynb | 111 ++++++++++++++++++++++++-- 1 file changed, 105 insertions(+), 6 deletions(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index 16fb4a718..e8786d0b6 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/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" ] }, { @@ -647,6 +648,7 @@ "# Compute divergence of the gradient (Laplacian)\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", @@ -759,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", @@ -834,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", @@ -907,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", From 9ab0388f077195161658da477f5f2579071fad55 Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Wed, 4 Nov 2026 09:59:55 -0500 Subject: [PATCH 7/8] fix small calculus .ipynb bugs, add 7y-10z example --- docs/user-guide/vector_calculus.ipynb | 63 +++++++++++++++++++++------ 1 file changed, 50 insertions(+), 13 deletions(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index e8786d0b6..923c30ffc 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/docs/user-guide/vector_calculus.ipynb @@ -811,7 +811,9 @@ "# 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", + "print(\n", + " f\"Max difference between laplacian and div_of_gradient: {np.abs(diff).max().item():.2e}\"\n", + ")\n", "laplacian" ] }, @@ -827,7 +829,7 @@ "$$\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." + "$x$, $y$, and $z$ can all be formed from linear combinations of spherical harmonics of degree 1, so taking the Laplacian of them, or any linear combination of them, should give results close to 0, aside from discretization errors." ] }, { @@ -840,22 +842,57 @@ "grid = uxds.uxgrid\n", "a = grid.sphere_radius\n", "\n", - "# Cartesian x of the face centers (unit sphere), wrapped as a UxDataArray\n", + "# Cartesian x, y, z 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", + "x_uarr = grid.face_x + zero\n", + "y_uarr = grid.face_y + zero\n", + "z_uarr = grid.face_z + zero\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", + "lap_x = x_uarr.laplacian()\n", + "residual_x = lap_x - (-2.0 * x_uarr / a**2)\n", + "finite_x = np.isfinite(residual_x) # exclude NaNs from comparisons below\n", + "res_finite_x = residual_x[finite_x]\n", + "\n", + "print(\n", + " f\" (x) Residual range: ({res_finite_x.min().item():.2e}, {res_finite_x.max().item():.2e})\"\n", + ")\n", + "# What does \"close to 0\" really mean? How close is \"close enough\"?\n", + "# Use max value of analytic solution as something to compare to:\n", + "max_analytic_x = (2 * np.abs(x_uarr[finite]).max() / a**2).item()\n", + "print(\n", + " f\" (x) Max residual is {res_finite_x.max().item() / max_analytic_x:.2%} \"\n", + " f\"of max |analytic| (={max_analytic_x:.2e})\"\n", + ")\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", + "# q = 7y - 10z is also a degree-1 combination of spherical harmonics, so laplacian(q) = -2q/a^2\n", + "q = 7 * y_uarr - 10 * z_uarr\n", + "lap_q = q.laplacian()\n", + "residual_q = lap_q - (-2.0 * q / a**2)\n", + "finite_q = np.isfinite(residual_q)\n", + "res_finite_q = residual_q[finite_q]\n", "\n", - "residual.plot(cmap=\"RdBu_r\", aspect=1).opts(\n", - " title=\"Laplacian of x minus (-2x/a²) (Should ≈ 0)\", colorbar=True\n", - ")" + "print(\n", + " f\"(7y-10z) Residual range: ({res_finite_q.min().item():.2e}, {res_finite_q.max().item():.2e})\"\n", + ")\n", + "max_analytic_q = (2 * np.abs(q.values[finite]).max() / a**2).item()\n", + "print(\n", + " f\"(7y-10z) Max residual is {res_finite_q.max().item() / max_analytic_q:.2%} \"\n", + " f\"of max |analytic| (={max_analytic_q:.2e})\"\n", + ")\n", + "\n", + "# plots of residuals in both cases:\n", + "plot_x = residual_x.plot(cmap=\"RdBu_r\", aspect=1).opts(\n", + " title=\"Laplacian of x minus (-2x/a²) (Should ≈ 0)\",\n", + " colorbar=True,\n", + " axiswise=True, # don't force same clim as plot_q\n", + ")\n", + "plot_q = residual_q.plot(cmap=\"RdBu_r\", aspect=1).opts(\n", + " title=\"Laplacian of q=(7y-10z) minus (-2q/a²) (Should ≈ 0)\",\n", + " colorbar=True,\n", + " axiswise=True, # don't force same clim as plot_x\n", + ")\n", + "plot_x + plot_q" ] }, { From a4bfd94706ca0573259b4a774d3b43aa81b3dd6f Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Wed, 4 Nov 2026 10:25:58 -0500 Subject: [PATCH 8/8] fix small typos (discovered mostly by copilot) --- docs/user-guide/vector_calculus.ipynb | 10 ++++++---- test/core/test_vector_calculus.py | 2 +- uxarray/core/dataarray.py | 8 ++++---- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/docs/user-guide/vector_calculus.ipynb b/docs/user-guide/vector_calculus.ipynb index 923c30ffc..0d714d3d6 100644 --- a/docs/user-guide/vector_calculus.ipynb +++ b/docs/user-guide/vector_calculus.ipynb @@ -824,9 +824,11 @@ "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", "$$\n", "\\nabla^2 q = -\\frac{l(l+1)}{a^2}\\, q\n", "$$\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, or any linear combination of them, should give results close to 0, aside from discretization errors." @@ -859,7 +861,7 @@ ")\n", "# What does \"close to 0\" really mean? How close is \"close enough\"?\n", "# Use max value of analytic solution as something to compare to:\n", - "max_analytic_x = (2 * np.abs(x_uarr[finite]).max() / a**2).item()\n", + "max_analytic_x = (2 * np.abs(x_uarr[finite_x]).max() / a**2).item()\n", "print(\n", " f\" (x) Max residual is {res_finite_x.max().item() / max_analytic_x:.2%} \"\n", " f\"of max |analytic| (={max_analytic_x:.2e})\"\n", @@ -875,7 +877,7 @@ "print(\n", " f\"(7y-10z) Residual range: ({res_finite_q.min().item():.2e}, {res_finite_q.max().item():.2e})\"\n", ")\n", - "max_analytic_q = (2 * np.abs(q.values[finite]).max() / a**2).item()\n", + "max_analytic_q = (2 * np.abs(q.values[finite_q]).max() / a**2).item()\n", "print(\n", " f\"(7y-10z) Max residual is {res_finite_q.max().item() / max_analytic_q:.2%} \"\n", " f\"of max |analytic| (={max_analytic_q:.2e})\"\n", @@ -1191,7 +1193,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -1205,7 +1207,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.2" + "version": "3.11.15" } }, "nbformat": 4, diff --git a/test/core/test_vector_calculus.py b/test/core/test_vector_calculus.py index 9d9e05918..60603af56 100644 --- a/test/core/test_vector_calculus.py +++ b/test/core/test_vector_calculus.py @@ -1115,7 +1115,7 @@ class TestLaplacianSphericalHarmonics: where q is any linear combination of any spherical harmonics of degree l. Linear combinations of harmonics of degree 1 include: x, y, z - Linear combinations of harmonics of degree 1 include: xy, xz, yz, x^2-y^2, 3z^2-r^2 + Linear combinations of harmonics of degree 2 include: xy, xz, yz, x^2-y^2, 3z^2-r^2 """ def grid_and_xyz(self, zoom): """returns (healpix grid, x, y, z) for this zoom. diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index 891a5c951..223dbf78a 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -1649,7 +1649,7 @@ def gradient(self, scale_by_radius: bool = True) -> UxDataset: the result carries physical units: [u units]/meter (e.g. result has units of 1/s if u represents velocity). When ``False`` (or if the grid's ``sphere_radius`` is set to 1.0), - the result caries per-radian units: [u units]/radian. + the result carries per-radian units: [u units]/radian. (e.g. result has units of (m/s)/rad if u represents velocity). CAUTION: this scaling does not account for any radius increases @@ -1839,7 +1839,7 @@ def divergence( the result carries physical units: [u units]/meter (e.g. result has units of 1/s if u represents velocity). When ``False`` (or if the grid's ``sphere_radius`` is set to 1.0), - the result caries per-radian units: [u units]/radian. + the result carries per-radian units: [u units]/radian. (e.g. result has units of (m/s)/rad if u represents velocity). CAUTION: this scaling does not account for any radius increases @@ -1961,9 +1961,9 @@ def laplacian(self, *, scale_by_radius: bool = True) -> "UxDataArray": Whether to divide unit-sphere derivatives by ``uxgrid.sphere_radius``. When ``True`` (and the grid's ``sphere_radius`` is set correctly), the result carries physical units: [q units]/meter^2 - (e.g. result has units of 1/s^2 if q represents velocity). + (e.g. result has units of 1/(m*s) if q represents velocity). When ``False`` (or if the grid's ``sphere_radius`` is set to 1.0), - the result caries per-radian^2 units: [q units]/radian^2. + the result carries per-radian^2 units: [q units]/radian^2. (e.g. result has units of (m/s)/rad^2 if q represents velocity). CAUTION: this scaling does not account for any radius increases