Skip to content

Clean up and vectorize the face_areas API (#1571)#1577

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

Clean up and vectorize the face_areas API (#1571)#1577
rajeeja wants to merge 3 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_dataarray flags.
  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 — 6.3x faster with identical areas.

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 [736bc03] After [84af9e9] Ratio Benchmark (Parameter)
- 580M 391M 0.67 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
- 699M 390M 0.56 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
- 499M 384M 0.77 mpas_ocean.Gradient.peakmem_gradient('480km')

Benchmarks that have stayed the same:

Change Before [736bc03] After [84af9e9] Ratio Benchmark (Parameter)
9.86±0.04μs 10.0±0.1μs 1.01 bench_connectivity.Connectivity.time_edge_face('120km')
10.4±0.2μs 10.3±0.1μs 0.99 bench_connectivity.Connectivity.time_edge_face('480km')
10.1±0.1μs 9.93±0.05μs 0.98 bench_connectivity.Connectivity.time_edge_node('120km')
10.1±0.08μs 10.3±0.04μs 1.02 bench_connectivity.Connectivity.time_edge_node('480km')
10.0±0.04μs 9.83±0.07μs 0.98 bench_connectivity.Connectivity.time_face_edge('120km')
10.1±0.08μs 10.4±0.04μs 1.03 bench_connectivity.Connectivity.time_face_edge('480km')
9.88±0.08μs 10.1±0.1μs 1.02 bench_connectivity.Connectivity.time_face_face('120km')
10.2±0.04μs 10.4±0.1μs 1.02 bench_connectivity.Connectivity.time_face_face('480km')
20.0±0.1μs 20.2±0.07μs 1.01 bench_connectivity.Connectivity.time_face_node('120km')
20.7±0.3μs 21.1±0.6μs 1.02 bench_connectivity.Connectivity.time_face_node('480km')
10.0±0.1μs 10.0±0.3μs 1 bench_connectivity.Connectivity.time_node_edge('120km')
10.1±0.1μs 10.2±0.1μs 1.01 bench_connectivity.Connectivity.time_node_edge('480km')
9.87±0.05μs 10.4±0.07μs 1.05 bench_connectivity.Connectivity.time_node_face('120km')
10.2±0.1μs 10.2±0.05μs 1 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'))
421M 421M 1 face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
13.0±0.1ms 12.9±0.06ms 0.99 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
3.42±0.09ms 3.42±0.09ms 1 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
17.5±0.1ms 17.4±0.09ms 0.99 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
2.02±0.05ms 2.04±0.04ms 1.01 face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
988±4ms 998±10ms 1.01 import.Imports.timeraw_import_uxarray
889±4ns 881±6ns 0.99 mpas_ocean.CheckNorm.time_check_norm('120km')
849±10ns 852±8ns 1 mpas_ocean.CheckNorm.time_check_norm('480km')
721±6ms 715±7ms 0.99 mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')
44.6±0.3ms 44.8±0.2ms 1.01 mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')
621±10μs 621±9μs 1 mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km')
556±10μs 547±8μs 0.98 mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')
5.10±0.04ms 5.05±0.02ms 0.99 mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('120km')
3.55±0.01ms 3.55±0.03ms 1 mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('480km')
3.37±0.01s 3.31±0.01s 0.98 mpas_ocean.ConstructFaceLatLon.time_welzl('120km')
212±0.9ms 214±0.6ms 1.01 mpas_ocean.ConstructFaceLatLon.time_welzl('480km')
19.7±0.01ms 19.7±0.02ms 1 mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')
1.13±0.04ms 1.12±0.02ms 0.99 mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')
10.6±0.02ms 10.6±0.02ms 1 mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')
732±20μs 752±10μs 1.03 mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')
673±20ms 653±4ms 0.97 mpas_ocean.CrossSections.time_const_lat('120km', 1)
338±5ms 331±2ms 0.98 mpas_ocean.CrossSections.time_const_lat('120km', 2)
173±4ms 172±1ms 0.99 mpas_ocean.CrossSections.time_const_lat('120km', 4)
489±2ms 488±2ms 1 mpas_ocean.CrossSections.time_const_lat('480km', 1)
247±1ms 245±1ms 0.99 mpas_ocean.CrossSections.time_const_lat('480km', 2)
127±2ms 128±1ms 1.01 mpas_ocean.CrossSections.time_const_lat('480km', 4)
22.0±0.03ms 22.0±0.05ms 1 mpas_ocean.DualMesh.time_dual_mesh_construction('120km')
2.55±0.02ms 2.57±0.02ms 1 mpas_ocean.DualMesh.time_dual_mesh_construction('480km')
833±3ms 843±5ms 1.01 mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)
49.3±0.4ms 49.8±1ms 1.01 mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)
72.4±0.2ms 74.4±0.6ms 1.03 mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)
5.44±0.1ms 5.35±0.05ms 0.98 mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)
404M 404M 1 mpas_ocean.Gradient.peakmem_gradient('120km')
164±0.3ms 165±0.5ms 1 mpas_ocean.Gradient.time_gradient('120km')
11.4±0.07ms 11.5±0.03ms 1.01 mpas_ocean.Gradient.time_gradient('480km')
192±2μs 189±1μs 0.98 mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('120km')
86.7±0.5μs 89.0±0.8μs 1.03 mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('480km')
351M 351M 1 mpas_ocean.Integrate.peakmem_integrate('120km')
330M 330M 1 mpas_ocean.Integrate.peakmem_integrate('480km')
182±2μs 178±1μs 0.98 mpas_ocean.Integrate.time_integrate('120km')
159±0.8μs 157±0.7μs 0.99 mpas_ocean.Integrate.time_integrate('480km')
188±1ms 189±1ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude')
188±0.7ms 188±1ms 1 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include')
188±1ms 189±2ms 1.01 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')
13.8±0.05ms 13.9±0.2ms 1.01 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude')
13.7±0.06ms 14.0±0.3ms 1.02 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include')
13.8±0.07ms 14.1±0.2ms 1.02 mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')
284±2μs 284±2μs 1 mpas_ocean.PointInPolygon.time_face_search_lonlat('120km')
285±1μs 285±1μs 1 mpas_ocean.PointInPolygon.time_face_search_lonlat('480km')
270±1μs 269±1μs 1 mpas_ocean.PointInPolygon.time_face_search_xyz('120km')
272±0.6μs 269±1μs 0.99 mpas_ocean.PointInPolygon.time_face_search_xyz('480km')
213±1ms 211±0.4ms 0.99 mpas_ocean.RemapDownsample.time_bilinear_remapping
226±0.9ms 228±0.6ms 1.01 mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping
4.26±0.04ms 4.26±0.01ms 1 mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping
1.18±0.01s 1.20±0.01s 1.01 mpas_ocean.RemapUpsample.time_bilinear_remapping
35.5±0.4ms 36.5±0.4ms 1.03 mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping
8.87±0.2ms 8.93±0.1ms 1.01 mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping
28.4±0.2ms 29.3±0.2ms 1.03 mpas_ocean.ZonalAverage.time_zonal_average('120km')
5.96±0.02ms 5.96±0.05ms 1 mpas_ocean.ZonalAverage.time_zonal_average('480km')
326M 326M 1 quad_hexagon.QuadHexagon.peakmem_open_dataset
325M 325M 1 quad_hexagon.QuadHexagon.peakmem_open_grid
6.46±0.1ms 6.46±0.2ms 1 quad_hexagon.QuadHexagon.time_open_dataset
5.48±0.07ms 5.48±0.1ms 1 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.

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

2 participants