warn if calculate_total_face_area() too large; add FESOM areas regression test

[Sevans711]

Closes #425

Overview

As discussed in thread in #425, the issue initially reported there has since been fixed by other changes. This PR adds a regression test (ensuring the total area for this FESOM grid is close to ~8.38 sr).

It also adds a warning to calculate_total_face_areas() if the result is larger than the total area of the sphere, to avoid silently returning wrong answers on malformed grids. It also updates the docstring for that function to clarify the warning might be raised.

Finally, a slight addition beyond just solving 425: updates the calculate_total_face_areas() docstring to clarify the actual behavior: it is not actually always equivalent to compute_face_areas().sum(). Specifically, for HEALPix grids, when using default values for all arguments, it is equivalent to face_areas.sum() (which respects HEALPix equal areas) instead.

Did not add any CI tests for this warning… please let me know if you think a regression test is necessary for it. Tested locally by temporarily locally adding result = result * 1.5 near the end of calculate_total_face_areas(), to ensure bad results actually lead to raising warning and failing FESOM area regression test.

PR Checklist

General

• An issue is created and linked
• Added appropriate labels (if your uxarray repo permissions allow it)
• Filled out Overview and Expected Usage (if applicable) sections

Testing & Benchmarking

• [N/A] Adequate tests are created if there is new functionality
• Tests are not too basic (such as simply calling a function and nothing else)
• Tests cover all major paths in your new functions
• [N/A] If this PR could affect performance, ran ASV benchmarks and confirmed they show expected behavior (add a new benchmark if necessary)

Documentation

• [N/A] Docstrings have been added to all new functions
• Docstrings have been updated with any function changes
• [N/A] User (public) functions have been added to docs/api.rst
• [N/A] Internal (private) function names start with an underscore (_)

[cmdupuis3]

Wouldn't it make sense to also check if result < ... * (1 - RTOL)?

[Sevans711]

That check would be like asking "is the total area less than the sphere's surface area?", right?

I wouldn't want a check like that because grids are always free to only cover a smaller region or have holes (e.g., see example in discussion of originally-linked issue).

[cmdupuis3]

Yeah. What I mean is that depending on how you approximate the sphere, you may be able to say something about the bounds, so it could be the case that this approximation is strictly larger than the sphere's total area. But if not, we should have tolerance on the lower bound also. But yeah, if we're not worried about preserving area on the lower bound, then we don't need to worry about it.

[Sevans711]

I think we aren't worried about the lower bound here, I'm mostly understanding the upper bound warning as a quick way to flag the possibility of a malformed grid.

Yeah. What I mean is that depending on how you approximate the sphere, you may be able to say something about the bounds, so it could be the case that this approximation is strictly larger than the sphere's total area. But if not, we should have tolerance on the lower bound also.

I'm not sure what this means, could you clarify a bit further (if still relevant, assuming you agree it is okay to not check a lower bound)?

[cmdupuis3]

So like, consider if you were trying to build a grid on the inside surface of the sphere, by placing nodes on the sphere. The area would approximate the sphere's area, but it would inherently be slightly smaller because it's strictly inside. In that case you'd expect the approximate area to be within (1 - RTOL) * sphere area. If we put faces tangent to the sphere, you'd have the opposite situation because the nodes are inherently outside of the sphere. What I'm getting at is depending on the method, you could have some mix of the two, or maybe you have some other way of knowing if the grid area is supposed to be strictly larger or smaller than the sphere.

[Sevans711]

Ahh, thank you, that helped me understand. I think that could be good to consider in the future, but isn't as relevant here. Here, maybe a clearer way to phrase the intention of the check would be to write something like:

full_sphere_area = 4 * np.pi * self.sphere_radius**2
assert (result < full_sphere_area) or np.isclose(result, full_sphere_area, rtol=RTOL)

I.e., the result is definitely small enough (less than full sphere area) or maybe this is a grid that covers the whole sphere (but account for possible rounding errors in calculations). Right now I am planning to keep the underlying code unchanged, but please let me know if you would want to see it rewritten or maybe with an extra comment added to clarify. (Note, the docstring does already include: Additionally, raises a warning if the result is larger than the total area of a sphere (4 * pi * self.sphere_radius**2).)