Compute `face_node_angles` from uxarray.Grid?

[Sevans711]

Proposed new feature or change:

This would have the same shape as face_node_connectivity: for every face, for every node in that face, this tells the angle at that node.

Computing angles is necessary for computing one common measure of grid skewness. Related: #1305.

There are a few different possible ways to get the result, so I'm not sure if this should be a cached property or a function with a few optional flags:

• should the result be in radians or degrees?
• should the computation use 'euclidean' / flat geometry (sum of angles of a triangle is always 180 degrees) or respect the 'spherical' geometry (sum of angles of a triangle is always larger than 180 degrees)?

What are the best defaults for those flags?

Also, perhaps grid.face_node_angles(geometry=...) should dispatch to face_node_angles_flat or face_node_angles_spherical properties (whose results should be cached)?

[Sevans711]

Pinging @rajeeja in case you have thoughts or advice on this!

[rajeeja]

Hi @Sevans711 — looked into this. Answers to your questions first, then a numerics point and a question for @hongyuchen1030.

Your questions

• Radians vs degrees: default radians (matches uxarray's internal convention across arcs.py / coordinates.py), with a degrees=True flag for convenience.
• Flat vs spherical: default spherical — it is the true geometry of these grids and the one that matters for skewness; keep euclidean as a cheaper option.
• Best defaults: geometry="spherical", degrees=False.
• Property vs function: go with grid.face_node_angles(geometry=..., degrees=...) dispatching to cached _face_node_angles_spherical / _face_node_angles_euclidean internals. Your instinct is right — a bare property cannot take flags, and the dispatch-to-cached-variant pattern gives both flexibility and caching.

Compute lazily, never eagerly. These angles are a derived, optional quantity most workflows never need, and computing them is a per-face/per-node pass over the whole mesh. Follow the face_areas pattern: compute on first request, cache, and never touch it at grid construction. Cache per geometry variant; keep the cache in radians and convert to degrees on return (cheap post-scale, no need to cache separately).

Geometry, so you don't compute the wrong thing. The angle at a node is the interior angle between the two edges meeting there, not the angle between node position vectors. At node B with face-neighbors A and C: euclidean is the angle between chords (A−B) and (C−B); spherical is the angle between the great-circle arcs BA and BC, from the tangent directions at B (project A and C onto the tangent plane at B, i.e. remove the component along B, then take the angle between the tangents) — equivalently the angle between the edge-plane normals n1 = B×A and n2 = B×C.

Reuse existing primitives: uxarray/grid/utils.py::_small_angle_of_2_vectors uses the stable 2*atan2(|u−v|,|u+v|) form — prefer it over anything arccos-based, since skewed / near-degenerate cells (your regime) are where arccos loses precision. node_x/y/z, face_node_connectivity, and n_nodes_per_face are cached and ready; mask INT_FILL_VALUE in the connectivity for variable-sided faces or the padding slots produce junk angles.

Free validation: for a spherical triangle, sum(interior angles) − π equals the spherical excess, which equals the face area — so you can cross-check against grid.face_areas in tests (generalizes to (n−2)π + area for n-gons; flat sum is exactly (n−2)π).

On accuracy / #1513. The spherical angle is built from cross products of the edge vectors (n1 = B×A, n2 = B×C). Naive cross products are difference-of-products and lose precision to cancellation exactly when edges are short or nearly collinear — which is the skewed-cell regime skewness analysis cares about most. The compensated cross-product kernels in the AccuSphGeom port (#1513) are built for precisely that, so I would structure the cross-product step so it can use those kernels once #1513 merges, rather than blocking on it now. #1513 is still open (not yet merged).

@hongyuchen1030 — a question for you, referencing this issue (#1566): AccuSphGeom's EFT is aimed at robust predicates (sign / orientation). Computing a node angle is a magnitude problem — accurate normals via accucross, then atan2(|n1×n2|, n1·n2). Does the compensated cross product meaningfully improve the accuracy of the resulting angle magnitude in the near-degenerate regime, or is the final error there dominated by the atan2 conditioning near 0/π? Trying to gauge whether wiring the #1513 kernels into this buys real angle accuracy or mainly helps the orientation/sign side.

Happy to review once you have a draft.

[Sevans711]

@rajeeja Thank you for your detailed response! My thoughts below:

Best defaults: geometry="spherical", degrees=False

Sounds good.

go with grid.face_node_angles(geometry=..., degrees=...) dispatching to cached _face_node_angles_spherical / _face_node_angles_euclidean internals.

Compute lazily, never eagerly. … keep the cache in radians.

Agreed.

Geometry [misc notes]

I think that all makes sense. I might realize I have questions once I start the implementation; I'll follow up if I get stuck.

Reuse existing primitives: uxarray/grid/utils.py::_small_angle_of_2_vectors

Thank you for pointing me to this! This helped me discover _angle_of_2_vectors, which is probably what I will end up using, to account for the possibility of non-convex faces where an internal angle can be >180 degrees. Unless… do uxarray grids actually require all faces to be convex?

Free validation: cross-check against grid.face_areas

I agree (and also, Claude happened to mention that idea as well when I initially asked about computing skewness). I found #1567 and #1571 while looking into face_areas because of this.

On accuracy: … [precision issues in the regime we care about, might be fixed but AccuSphGeom port]

This makes sense, and I will keep in mind that we might replace with more accurate cross products, or maybe a flag to swap between modes, in the near future.

Happy to review once you have a draft.

Thank you! I will request you as a reviewer once the PR is ready to go.

[erogluorhan]

Some thoughts:

• Computation of this should be straight forward and inexpensive since we already have face lats/lons and node lats/lons
• Calculations should be on spherical geometry only I believe
• Caching would be costly and redundant though - I think we should avoid it
  • We had some initial discussions about similar things and ended up deciding to cache only the necessary-to-store properties and deferring to adhoc calculations for other properties like this. That approach still makes sense to me from a memory utilization standpoint.

[Sevans711]

@erogluorhan thank you for sharing these thoughts as well. I'll focus on spherical geometry first, and disable caching by default. I might still add a flag like grid.face_node_angles(..., cache=False) to enable the user to opt-in to caching, in case they know ahead of time that they will be doing multiple computations related to angles, if that sounds reasonable?

[erogluorhan]

@erogluorhan thank you for sharing these thoughts as well. I'll focus on spherical geometry first

Great; thank you!

, and disable caching by default. I might still add a flag like grid.face_node_angles(..., cache=False) to enable the user to opt-in to caching, in case they know ahead of time that they will be doing multiple computations related to angles, if that sounds reasonable?

This is worth running some large data tests to see the performance and memory usage on NCAR's HPC clusters, maybe using 3.75 km data. Memory usage associated with the caching of it could be roughly calculated actually. Assuming you'll use float, it would (for our sample 3.75 km MPAS grid) be around (~42M faces * 6 nodes * 8 bytes) = ~2GB of extra memory usage.

The rest should be:

• If execution time is negligibly short (I am anticipating this to be the case), we should definitely avoid caching
• Else, let's further discuss expectations vs costs with caching

[Sevans711]

Small update: I managed to finish an initial implementation (see commit ecba4f1). Still need to add tests and examples before opening a PR.

The non-spherical geometry seems simple to support (it doesn't add any operations, it just skips the "project onto the tangent plane" step) so I included it as an option for now; it only took a few extra minutes to implement, to document the relevant kwarg and add "if geometry=..." checks inside the code.

I definitely agree that caching should be disabled by default, unless the calculation actually becomes a bottleneck. I did still include the caching option though, because it only took a few extra minutes to add it. This should help us later test in real workflows whether it saves non-negligible amounts of time:

1. Run workflow/analysis without ever touching the cache kwarg. This will avoid caching entirely.
2. Run same workflow/analysis but call grid.face_node_angles(..., cache=True) once at the top. This will lead to all calls using the cached value.

If tests show there's never any point in supporting caching and we want to remove that option to avoid the overhead of slightly longer code / more for users to read about and consider, it should be easy enough to remove the relevant pieces of the code.

[erogluorhan]

Re caching: Please see my previous comment and complete the tests accordingly and let us know the results before introducing functionality into API. (Once introduced into user API, even if it is disabled by default, removing it later would break user workflows if it gets some usage)

Re flat vs sphere: @rajeeja @hongyuchen1030 do you think it is worth to introduce any flat geometry calculations or should we focus only on spherical geometry?

[erogluorhan]

Small update: I managed to finish an initial implementation (see commit ecba4f1). Still need to add tests and examples before opening a PR.

Oh, I just had a look into the code now, and I believe it could be and has to be vectorized (instead of nested for loops) to run it in virtually "no time" for varying grid sizes. That way, I believe we'd have negligible execution times with it.