From 2ab1489f84e21e0a4c3ba26c9c9d6d909f24a073 Mon Sep 17 00:00:00 2001
From: Rajeev Jain <rajeeja@gmail.com>
Date: Thu, 6 Aug 2026 09:34:27 -0500
Subject: [PATCH 1/2] Restore cached face_areas fast path in
 calculate_total_face_area

The HEALPix fast path assigned result but did not return, so it was
always overwritten by the recomputed quadrature sum. Restores the
else branch and tightens RTOL back to 1e-7.
---
 test/grid/grid/test_areas.py | 27 +++++++++++++++++++++++++++
 uxarray/grid/grid.py         | 16 ++++++++--------
 2 files changed, 35 insertions(+), 8 deletions(-)

diff --git a/test/grid/grid/test_areas.py b/test/grid/grid/test_areas.py
index f11cbde8d..f1d99f74f 100644
--- a/test/grid/grid/test_areas.py
+++ b/test/grid/grid/test_areas.py
@@ -169,3 +169,30 @@ def test_face_areas_fesom(gridpath):
     uxgrid = ux.open_grid(gridpath("ugrid", "fesom", "fesom.mesh.diag.nc"))
     total_area = uxgrid.calculate_total_face_area()
     nt.assert_almost_equal(total_area, 8.3780, decimal=4)
+
+
+def test_total_face_area_healpix_uses_cached_equal_areas():
+    """Default args must reuse the cached ``face_areas``, not recompute.
+
+    HEALPix faces are exactly equal-area, so the cached sum is exactly 4*pi.
+    Recomputing via quadrature drifts and loses that property.
+    """
+    uxgrid = ux.Grid.from_healpix(zoom=2)
+
+    total_area = uxgrid.calculate_total_face_area()
+
+    nt.assert_allclose(total_area, np.sum(uxgrid.face_areas.values), rtol=0)
+    nt.assert_allclose(total_area, 4 * np.pi, rtol=1e-12)
+
+
+def test_total_face_area_honors_quadrature_kwargs():
+    """Non-default quadrature settings must still trigger a fresh computation."""
+    uxgrid = ux.Grid.from_healpix(zoom=2)
+
+    recomputed = uxgrid.calculate_total_face_area(quadrature_rule="gaussian", order=2)
+
+    nt.assert_allclose(
+        recomputed,
+        np.sum(uxgrid.compute_face_areas(quadrature_rule="gaussian", order=2)),
+        rtol=0,
+    )
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 779367743..67c178c7a 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -2015,16 +2015,16 @@ def calculate_total_face_area(
             and not latitude_adjusted_area
         ):
             result = np.sum(self.face_areas.values)
-
-        result = np.sum(
-            self.compute_face_areas(
-                quadrature_rule=quadrature_rule,
-                order=order,
-                latitude_adjusted_area=latitude_adjusted_area,
+        else:
+            result = np.sum(
+                self.compute_face_areas(
+                    quadrature_rule=quadrature_rule,
+                    order=order,
+                    latitude_adjusted_area=latitude_adjusted_area,
+                )
             )
-        )
 
-        RTOL = 1e-6  # 1e-7 had warnings in existing CI tests (as of 2026-08-05), 1e-6 did not.
+        RTOL = 1e-7
         if result > 4 * np.pi * self.sphere_radius**2 * (1 + RTOL):
             warnings.warn(
                 f"Total face area (={result}) exceeds the surface area of the whole sphere "

From 4356eba2198e4d00303b8bb36acaa51207035c95 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Thu, 6 Aug 2026 15:33:53 -0400
Subject: [PATCH 2/2] comment to say RTOL face areas is mostly arbitrary

---
 uxarray/grid/grid.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 67c178c7a..2f1bb582d 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -2024,6 +2024,8 @@ def calculate_total_face_area(
                 )
             )
 
+        # Choose RTOL. Mostly just an arbitrary decision....
+        # but noting that 1e-9 had warnings in existing CI tests (as of 2026-08-06), while 1e-8 did not.
         RTOL = 1e-7
         if result > 4 * np.pi * self.sphere_radius**2 * (1 + RTOL):
             warnings.warn(