Skip to content

Clean up and vectorize the face_areas API#1577

Open
rajeeja wants to merge 8 commits into
mainfrom
rajeeja/area_cleanup
Open

Clean up and vectorize the face_areas API#1577
rajeeja wants to merge 8 commits into
mainfrom
rajeeja/area_cleanup

Conversation

@rajeeja

@rajeeja rajeeja commented Jul 15, 2026

Copy link
Copy Markdown
Contributor

Closes #1571.

  1. Un-deprecate compute_face_areas() as the public quadrature-capable entry point — returns areas by default, with return_jacobian and as_uxarray flags (the latter returns a UxDataArray paired with the grid).
  2. Rename the internal worker to _compute_face_areas_and_jacobian and privatize the low-level area.py helpers.
  3. Document the quadrature_rule/order/latitude_adjusted_area options; note calculate_total_face_area is equivalent to compute_face_areas().sum().
  4. Vectorize get_all_face_area_from_coords with prange and hoist the per-face quadrature setup out of the loop — ~6x faster with identical areas; adds a FaceAreas asv benchmark.
  5. Plot face areas in area_calc.ipynb and the geometric/theoretical ratio in healpix.ipynb using the new as_uxarray flag.

Un-deprecate compute_face_areas() as the public quadrature-capable entry
point (areas by default, return_jacobian and as_dataarray flags), rename
the internal worker to _compute_face_areas_and_jacobian, privatize the
low-level area.py helpers, and document the quadrature options. Vectorize
get_all_face_area_from_coords with prange and hoist the per-face
quadrature setup out of the loop for a 6.3x speedup with identical areas.
Closes #1571.
@review-notebook-app

Copy link
Copy Markdown

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@Sevans711 Sevans711 added documentation Improvements or additions to documentation redesign Content relating to the ongoing redesign improvement Improvements on existing features or infrastructure run-benchmark Run ASV benchmark workflow labels Jul 15, 2026
@Sevans711
Sevans711 self-requested a review July 15, 2026 23:16
@Sevans711

Copy link
Copy Markdown
Collaborator

Thank you for these changes! I added run-benchmarks label but then realized I'm not sure if there are any benchmarks which actually rely on face_areas. If there's no visible impact, then would you be able to add at least one face_areas benchmark? If the impact is visible already, I don't feel strongly about needing a face_areas specific benchmark.

I'll come back to this soon (probably tomorrow), and if the benchmark impact is visible then I will review the rest!

@github-actions

github-actions Bot commented Jul 15, 2026

Copy link
Copy Markdown

ASV Benchmarking

Benchmark Comparison Results

Benchmarks that have improved:

Change Before [b0487c4] After [5025ddb] Ratio Benchmark (Parameter)
- 579M 391M 0.67 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
- 701M 390M 0.56 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
- 487M 384M 0.79 mpas_ocean.FaceAreas.peakmem_compute_face_areas('480km')
- 157±0.4ms 57.3±0.09ms 0.36 mpas_ocean.FaceAreas.time_compute_face_areas('120km')
- 12.0±0.08ms 5.93±0.05ms 0.49 mpas_ocean.FaceAreas.time_compute_face_areas('480km')
- 490M 384M 0.78 mpas_ocean.Gradient.peakmem_gradient('480km')

Benchmarks that have stayed the same:

Change Before [b0487c4] After [5025ddb] Ratio Benchmark (Parameter)
9.98±0.09μs 10.1±0.1μs 1.01 bench_connectivity.Connectivity.time_edge_face('120km')
10.4±0.2μs 10.1±0.05μs 0.98 bench_connectivity.Connectivity.time_edge_face('480km')
9.94±0.09μs 9.92±0.02μs 1 bench_connectivity.Connectivity.time_edge_node('120km')
10.1±0.1μs 10.1±0.08μs 1 bench_connectivity.Connectivity.time_edge_node('480km')
9.93±0.07μs 10.0±0.03μs 1.01 bench_connectivity.Connectivity.time_face_edge('120km')
10.4±0.2μs 10.2±0.1μs 0.98 bench_connectivity.Connectivity.time_face_edge('480km')
10.2±0.2μs 9.92±0.07μs 0.97 bench_connectivity.Connectivity.time_face_face('120km')
10.1±0.09μs 10.2±0.08μs 1.01 bench_connectivity.Connectivity.time_face_face('480km')
20.3±0.2μs 21.7±2μs 1.07 bench_connectivity.Connectivity.time_face_node('120km')
20.5±0.2μs 20.5±0.2μs 1 bench_connectivity.Connectivity.time_face_node('480km')
10.1±0.06μs 9.86±0.05μs 0.98 bench_connectivity.Connectivity.time_node_edge('120km')
10.4±0.3μs 10.1±0.05μs 0.97 bench_connectivity.Connectivity.time_node_edge('480km')
10.0±0.1μs 9.81±0.06μs 0.98 bench_connectivity.Connectivity.time_node_face('120km')
10.3±0.2μs 10.2±0.07μs 0.99 bench_connectivity.Connectivity.time_node_face('480km')
389M 389M 1 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
419M 419M 1 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
12.8±0.05ms 12.7±0.04ms 0.99 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
3.41±0.09ms 3.28±0.07ms 0.96 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
17.4±0.07ms 17.4±0.08ms 1 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
2.00±0.02ms 2.01±0.04ms 1 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
987±4ms 990±4ms 1 import.Imports.timeraw_import_uxarray
884±10ns 892±3ns 1.01 mpas_ocean.CheckNorm.time_check_norm('120km')
885±4ns 883±30ns 1 mpas_ocean.CheckNorm.time_check_norm('480km')
713±2ms 729±4ms 1.02 mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')
44.6±0.2ms 44.9±0.1ms 1.01 mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')
634±5μs 635±10μs 1 mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km')
550±10μs 549±4μs 1 mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')
5.07±0.04ms 5.05±0.03ms 1 mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('120km')
3.51±0.03ms 3.54±0.03ms 1.01 mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('480km')
3.30±0.03s 3.29±0.01s 1 mpas_ocean.ConstructFaceLatLon.time_welzl('120km')
218±2ms 215±4ms 0.99 mpas_ocean.ConstructFaceLatLon.time_welzl('480km')
19.7±0.02ms 19.7±0.02ms 1 mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')
1.12±0.01ms 1.11±0.03ms 1 mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')
10.6±0.03ms 10.6±0.02ms 1 mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')
741±20μs 731±20μs 0.99 mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')
657±3ms 651±0.7ms 0.99 mpas_ocean.CrossSections.time_const_lat('120km', 1)
331±0.4ms 329±0.4ms 1 mpas_ocean.CrossSections.time_const_lat('120km', 2)
170±0.7ms 171±0.5ms 1.01 mpas_ocean.CrossSections.time_const_lat('120km', 4)
489±2ms 487±3ms 1 mpas_ocean.CrossSections.time_const_lat('480km', 1)
245±0.6ms 247±1ms 1.01 mpas_ocean.CrossSections.time_const_lat('480km', 2)
126±0.9ms 128±0.4ms 1.01 mpas_ocean.CrossSections.time_const_lat('480km', 4)
22.0±0.05ms 22.0±0.03ms 1 mpas_ocean.DualMesh.time_dual_mesh_construction('120km')
2.56±0.01ms 2.59±0.02ms 1.01 mpas_ocean.DualMesh.time_dual_mesh_construction('480km')
399M 400M 1 mpas_ocean.FaceAreas.peakmem_compute_face_areas('120km')
839±7ms 838±7ms 1 mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)
49.5±0.7ms 50.4±0.6ms 1.02 mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)
72.5±0.2ms 72.3±0.5ms 1 mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)
5.45±0.04ms 5.37±0.05ms 0.99 mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)
404M 404M 1 mpas_ocean.Gradient.peakmem_gradient('120km')
164±0.4ms 165±0.5ms 1.01 mpas_ocean.Gradient.time_gradient('120km')
11.7±0.2ms 11.4±0.03ms 0.98 mpas_ocean.Gradient.time_gradient('480km')
191±1μs 190±2μs 0.99 mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('120km')
87.7±0.8μs 87.0±0.8μs 0.99 mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('480km')
351M 355M 1.01 mpas_ocean.Integrate.peakmem_integrate('120km')
330M 331M 1 mpas_ocean.Integrate.peakmem_integrate('480km')
181±0.4μs 178±3μs 0.98 mpas_ocean.Integrate.time_integrate('120km')
162±2μs 160±0.7μs 0.99 mpas_ocean.Integrate.time_integrate('480km')
189±3ms 186±2ms 0.98 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude')
186±1ms 187±0.9ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include')
187±1ms 185±1ms 0.99 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')
13.8±0.08ms 13.7±0.06ms 0.99 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude')
13.6±0.07ms 13.9±0.2ms 1.02 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include')
13.8±0.05ms 13.6±0.06ms 0.99 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')
286±1μs 285±2μs 1 mpas_ocean.PointInPolygon.time_face_search_lonlat('120km')
285±2μs 283±0.8μs 0.99 mpas_ocean.PointInPolygon.time_face_search_lonlat('480km')
270±1μs 272±2μs 1 mpas_ocean.PointInPolygon.time_face_search_xyz('120km')
270±1μs 270±3μs 1 mpas_ocean.PointInPolygon.time_face_search_xyz('480km')
212±0.2ms 213±0.8ms 1 mpas_ocean.RemapDownsample.time_bilinear_remapping
228±1ms 230±0.9ms 1.01 mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping
4.25±0.02ms 4.23±0.04ms 0.99 mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping
1.19±0s 1.18±0s 1 mpas_ocean.RemapUpsample.time_bilinear_remapping
35.4±0.5ms 36.3±0.1ms 1.03 mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping
9.03±0.2ms 8.79±0.03ms 0.97 mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping
29.0±0.2ms 29.4±0.4ms 1.02 mpas_ocean.ZonalAverage.time_zonal_average('120km')
6.10±0.05ms 5.94±0.02ms 0.97 mpas_ocean.ZonalAverage.time_zonal_average('480km')
326M 328M 1.01 quad_hexagon.QuadHexagon.peakmem_open_dataset
326M 326M 1 quad_hexagon.QuadHexagon.peakmem_open_grid
6.50±0.09ms 6.52±0.1ms 1 quad_hexagon.QuadHexagon.time_open_dataset
5.46±0.06ms 5.50±0.1ms 1.01 quad_hexagon.QuadHexagon.time_open_grid

@Sevans711 Sevans711 left a comment

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.

Looking through the plan of action discussed in #1571:

  1. (Looks good to me) Un-deprecate and make compute_face_areas() the public, kwargs-capable entry point. Return face_areas only by default, with return_jacobian=False (numpy-style) for callers that also want the jacobian.
  2. (Looks good to me) Rename _compute_face_areas_compute_face_areas_and_jacobian
  3. (Not completed as planned; see notes below) Add an optional as_dataarray=False flag to return a UxDataArray (result.uxgrid = self) so per-face areas are easy to plot.
  4. (Looks good to me) Keep calculate_total_face_area() but document it as equivalent to compute_face_areas().sum().
  5. (Looks good to me) Mark the low-level area.py routines private
  6. (Mostly good but I have more notes below) Vectorize the core computation.
  7. (Mostly good but I have more notes below) Fix the user guide (area_calc.ipynb).
  8. (Looks good, I did a small commit to help) improve docstrings for face_areas property and compute_face_areas().

More notes / requests / suggestions:

  • (3) Originally we discussed returning a UxDataArray object, not an xarray.DataArray. The current version here returns an xarray.DataArray. Request: make the optional flag return UxDataArray instead. I think UxDataArray is the way to go; it is more convenient, and I think the whole point for this optional flag is to return something more convenient.
    • For example, plotting result from current version looks something like: ux.UxDataArray(uxds.uxgrid.compute_face_areas(as_dataarray=True), uxgrid=uxds.uxgrid).plot(). But, with this requested change, it would become simpler: uxds.uxgrid.compute_face_areas(as_dataarray=True).plot().
    • I wonder if as_dataarray is the wrong name here... Request/suggestion: rename the flag to as_uxarray or as_uxdataarray to clarify it is a UxDataArray instead? Of the three names, my preferences would be as_uxarray (concise & clear) > as_uxdataarray (verbose & clear) > as_dataarray (sounds like xarray not uxarray). (Let me know if you want me to take on the task of setting this up!)
  • (7a) Request: plot face_areas somewhere in area_calc.ipynb, using the new optional flag. Either in a new section, or inside an existing section, e.g. maybe it makes sense to add in the "Calculate Area from Multiple Faces in Spherical Coordinates" section?
  • (7b) In the healpix.ipynb guide, there is a section comparing geometric face area calculations versus theoretical healpix areas. This has been updated to now be equivalent to grid.compute_face_faces() / grid.face_areas. There is also some numerical analysis of how the results vary across the grid. Request: plot this face_areas ratio across the grid. After making changes in response to point (3) above, it should be as simple as: (grid.compute_face_faces(as_uxarray=True) / grid.face_areas).plot().
  • (6) I also measured roughly 6x faster face_areas computations on main versus on this branch. This is a great improvement! Some of the ASV benchmarks show improvements in reducing peak memory usage, but none of them show this significant speedup. Request: add at least one benchmark that gets improved by these changes.

@erogluorhan erogluorhan changed the title Clean up and vectorize the face_areas API (#1571) Clean up and vectorize the face_areas API Jul 17, 2026
@Sevans711

Copy link
Copy Markdown
Collaborator

@rajeeja thank you for adding the total face area error check. I don't want to interrupt you if you're still working on things, but wanted to share my thoughts now while I am thinking about it:

  1. Maybe this check should only be in calculate_total_face_areas(), or in grid.validate()? It seems unusual to impose inefficiency on the low-level _compute_face_areas_and_jacobian method for a validation check.
  2. Flagging that face_areas will be much much larger than 4 * pi for any grids which are not on the unit sphere! From skimming your changes, it is unclear to me whether that would trigger this warning or not. Example non-unit-sphere-grid: 'uxarray/test/meshfiles/mpas/dyamond-30km/gradient_grid_subset.nc'

compute_face_areas() gains as_uxarray (renamed from as_dataarray) returning
a UxDataArray paired with the grid; plot face areas in area_calc.ipynb and
the geometric/theoretical ratio in healpix.ipynb using it; add a FaceAreas
asv benchmark capturing the vectorization speedup.
@rajeeja

rajeeja commented Jul 17, 2026

Copy link
Copy Markdown
Contributor Author

@Sevans711 thanks for the review — fixed the items: compute_face_areas() now returns a UxDataArray via as_uxarray=True (renamed from as_dataarray); face areas plotted in area_calc.ipynb and the geometric/theoretical ratio in healpix.ipynb using it; and added a FaceAreas asv benchmark that captures the speedup.

@rajeeja
rajeeja requested a review from Sevans711 July 17, 2026 18:54
@rajeeja

rajeeja commented Jul 17, 2026

Copy link
Copy Markdown
Contributor Author

@Sevans711 Good catches on both. You're right the check misfires on grids whose face_areas are parsed from the file in physical units — e.g. gradient_grid_subset.nc returns ~1.5e11 m² and would falsely warn, since the property returns file values, not unit-sphere steradians. There's no reliable way to know the unit at that level, so I've removed the check. #425's actual issue (the wrong value) is already fixed independently; I'll leave #425 open for a proper unit-aware validation later rather than ship a check that false-positives.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

documentation Improvements or additions to documentation improvement Improvements on existing features or infrastructure redesign Content relating to the ongoing redesign run-benchmark Run ASV benchmark workflow

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Clean up face_areas API

3 participants