-
Notifications
You must be signed in to change notification settings - Fork 45
Zonal mean helpers Fix #955
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Zonal mean helpers Fix #955
Conversation
amberchen122
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks good to me!
|
The
May you investigate the grids that are failing?
|
It looks like the zonal_mean branch also update the face-bounds function. I am looking into it right now. |
|
I’ve pinpointed where the issue occurs, and it’s once again related to the unnormalized coordinates in MPAS. Specifically, it's connected to the file In the def extreme_gca_latitude(gca_cart, extreme_type):
extreme_type = extreme_type.lower()
if extreme_type not in ("max", "min"):
raise ValueError("extreme_type must be either 'max' or 'min'")
n1, n2 = gca_cart
dot_n1_n2 = np.dot(n1, n2)
denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0)
d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom
d_a_max = (
np.clip(d_a_max, 0, 1)
if np.isclose(d_a_max, [0, 1], atol=MACHINE_EPSILON).any()
else d_a_max
)
# Before we make sure the grid coordinates are normalized, do not try to skip the normalization steps!
_, lat_n1 = _xyz_to_lonlat_rad_no_norm(n1[0], n1[1], n1[2])
_, lat_n2 = _xyz_to_lonlat_rad_no_norm(n2[0], n2[1], n2[2])
if 0 < d_a_max < 1:
node3 = (1 - d_a_max) * n1 + d_a_max * n2
node3 = np.array(_normalize_xyz_scalar(node3[0], node3[1], node3[2]))
d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1))
return (
max(d_lat_rad, lat_n1, lat_n2)
if extreme_type == "max"
else min(d_lat_rad, lat_n1, lat_n2)
)
else:
return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2)The latest changes are highly risky if we haven’t ensured that the input grid coordinates are normalized uniformly: _, lat_n1 = _xyz_to_lonlat_rad_no_norm(n1[0], n1[1], n1[2])
_, lat_n2 = _xyz_to_lonlat_rad_no_norm(n2[0], n2[1], n2[2])Because if n1 and n2 are not normalized, the returned longitude won’t be in the range of [0,360], and the latitude won’t be within [0,90] or [0,-90] degrees, which will lead to the corruption of all subsequent functions. This makes it critical to be extremely cautious when using |
|
Now, I have added the following test in gridfile_CSne8 = current_path / "meshfiles" / "scrip" / "outCSne8" / "outCSne8.nc"
grid_quad_hex = current_path/ "meshfiles" / "ugrid" / "quad-hexagon" / "grid.nc"
grid_geoflow = current_path/ "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
grid_mpas = current_path/ "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
grid_files_latlonBound = [grid_quad_hex, grid_geoflow, gridfile_CSne8, grid_mpas]
class TestLatlonBoundsFiles:
def test_face_bounds(self):
"""Test to ensure ``Grid.face_bounds`` works correctly for all grid
files."""
for grid_path in grid_files_latlonBound:
try:
# Open the grid file
self.uxgrid = ux.open_grid(grid_path)
# Test: Ensure the bounds are obtained
bounds = self.uxgrid.bounds
assert bounds is not None, f"Grid.face_bounds should not be None for {grid_path}"
except Exception as e:
# Print the failing grid file and re-raise the exception
print(f"Test failed for grid file: {grid_path}")
raise e
finally:
# Clean up the grid object
del self.uxgrid |
philipc2
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've made a change to asv.conf.json that should fix the failing tests (they were due to the following)
Co-authored-by: Philip Chmielowiec <[email protected]>
|
Benchmarks that have got worse:
I will profile this current implementation and see what we need to adjust tommorow. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I believe the root cause of the decrease in performance was that most of the existing calls to our numba wrapped numpy functions were removed.
Please replace any the following functions with their respective ones from uxarray.utils.computing
np.allnp.isclosenp.allclosenp.crossnp.dot
uxarray/uxarray/utils/computing.py
Lines 7 to 61 in c66286b
| @njit | |
| def all(a): | |
| """Numba decorated implementation of ``np.all()`` | |
| See Also | |
| -------- | |
| numpy.all | |
| """ | |
| return np.all(a) | |
| @njit | |
| def isclose(a, b, rtol=1e-05, atol=1e-08): | |
| """Numba decorated implementation of ``np.isclose()`` | |
| See Also | |
| -------- | |
| numpy.isclose | |
| """ | |
| return np.isclose(a, b, rtol=rtol, atol=atol) | |
| @njit | |
| def allclose(a, b, rtol=1e-05, atol=1e-08): | |
| """Numba decorated implementation of ``np.allclose()`` | |
| See Also | |
| -------- | |
| numpy.allclose | |
| """ | |
| return np.allclose(a, b, rtol=rtol, atol=atol) | |
| @njit | |
| def cross(a, b): | |
| """Numba decorated implementation of ``np.cross()`` | |
| See Also | |
| -------- | |
| numpy.cross | |
| """ | |
| return np.cross(a, b) | |
| @njit | |
| def dot(a, b): | |
| """Numba decorated implementation of ``np.dot()`` | |
| See Also | |
| -------- | |
| numpy.dot | |
| """ | |
| return np.dot(a, b) |
I have updated the current implementation with their respective ones from I did found some leftover calls in the main branch, but we can fix them after we merge the zonal mean |
|
Looks like there is still about a Benchmarks that have improved:
Benchmarks that have stayed the same:
Benchmarks that have got worse:
|
|
Latest benchmark runs. After a few changes, the performance is only slightly worse, which is acceptable and will be ironed out in Will make a few final changes and get this merged today. |
philipc2
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Nice work!
Thanks for your great work! |

Overview
main._get_zonal_faces_weight_at_constLatgca_constlat_intersection_pole_point_inside_polygonpoint_within_gcafunction used in the Zonal Mean PR.Expected Usage
PR Checklist
General
Testing
Documentation
_) and have been added todocs/internal_api/index.rstdocs/user_api/index.rstExamples
docs/examples/folderdocs/examples.rsttoctreedocs/gallery.ymlwith appropriate thumbnail photo indocs/_static/thumbnails/